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Foreword 


Dear reader, 

Our aim with the series Simula SpringerBriefs on Computing is to provide 
compact introductions to selected fields of computing. Entering a new field of 
research can be quite demanding for graduate students, postdocs, and experienced 
researchers alike: the process often involves reading hundreds of papers, and the 
methods, results and notation styles used often vary considerably, which makes for 
a time-consuming and potentially frustrating experience. The briefs in this series are 
meant to ease the process by introducing and explaining important concepts and 
theories in a relatively narrow field, and by posing critical questions on the fun- 
damentals of that field. A typical brief in this series should be around 100 pages and 
should be well suited as material for a research seminar in a well-defined and 
limited area of computing. 

We have decided to publish all items in this series under the SpringerOpen 
framework, as this will allow authors to use the series to publish an initial version 
of their manuscript that could subsequently evolve into a full-scale book on a 
broader theme. Since the briefs are freely available online, the authors will not 
receive any direct income from the sales; however, remuneration is provided for 
every completed manuscript. Briefs are written on the basis of an invitation from a 
member of the editorial board. Suggestions for possible topics are most welcome 
and can be sent to aslak @simula.no. 


January 2016 Prof. Aslak Tveito 
CEO 


Dr. Martin Peters 
Executive Editor Mathematics 
Springer Heidelberg, Germany 


Preface 


Finding proper values of physical parameters in mathematical models is often 
quite a challenge. While many have gotten away with using just the math- 
ematical symbols when doing science and engineering with pen and paper, 
the modern world of numerical computing requires each physical parameter 
to have a numerical value, otherwise one cannot get started with the com- 
putations. For example, in the simplest possible transient heat conduction 
simulation, a case relevant for a real physical material needs values for the 
heat capacity, the density, and the heat conduction coefficient of the ma- 
terial. In addition, relevant values must be chosen for initial and boundary 
temperatures as well as the size of the material. With a dimensionless math- 
ematical model, as explained in Chapter 3.2, no physical quantities need to 
be assigned (!). Not only is this a simplification of great convenience, as one 
simulation is valid for any type of material, but it also actually increases the 
understanding of the physical problem. 

Scaling of differential equations is basically a simple mathematical process, 
consisting of the chain rule for differentiation and some algebra. The choice 
of scales, however, is a non-trivial topic, which may cause confusion among 
practitioners without extensive experience with scaling. How to choose scales 
is unfortunately not well treated in the literature. Most of the times, authors 
just state scales without proper motivation. The choice of scales is highly 
problem-dependent and requires knowledge of the characteristic features of 
the solution or the physics of the problem. The present notes aim at explaining 
“all nuts and bolts” of the scaling technique, including choice of scales, the 
algebra, the interpretation of dimensionless parameters in scaled models, and 
how scaling impacts software for solving differential equations. 

Traditionally, scaling was mainly used to identify small parameters in 
mathematical models, such that perturbation methods based on series ex- 
pansions in terms of the small parameters could be used as an approximate 
solution method for differential equations. Nowadays, the greatest practical 
benefit of scaling is related to running numerical simulations, since scaling 
greatly simplifies the choice of values for the input data and makes the sim- 
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ulations results more widely applicable. The number of parameters in scaled 
models may be much less than the number of physical parameters in the 
original model. The parameters in scaled models are also dimensionless and 
express ratios of physical effects rather than levels of individual effects. Set- 
ting meaningful values of a few dimensionless numbers is much easier than 
determining physically relevant values for the original physical parameters. 

Another great benefit of scaling is the physical insight that follows from 
dimensionless parameters. Since physical effects enter the problem through a 
few dimensionless groups, one can from these groups see how different effects 
compete in their impact on the solution. Ideally, a good physical understand- 
ing should provide the same insight, but it is not always easy to “think right” 
and realize how spatial and temporal scales interact with physical parame- 
ters. This interaction becomes clear through the dimensionless numbers, and 
such numbers are therefore a great help, especially for students, in developing 
a correct physical understanding. 

Since we have a special focus on scaling related to numerical simulations, 
the notes contain a lot of examples on how to program with dimensionless 
differential equation models. Most numerical models feature quantities with 
dimension, so we show in particular how to utilize such existing models to 
solve the equations in the associated scaled model. 

Scaling is not a universal mathematical technique as the details depend 
on the problem at hand. We therefore present scaling in a range of specific 
applications, starting with simple ODEs, progressing with basic PDEs, before 
attacking more complicated models, especially from fluid mechanics. 

Chapter 1 discusses units and how to make programs that can automat- 
ically take care of unit conversion (the most frequent mathematical mistake 
in industry and science?). Section 2.1 introduces the mathematics of scaling 
and the thinking about scales in a simple ODE problem modeling expo- 
nential decay. The ideas are generalized to nonlinear ODEs and to systems 
of ODEs. Another ODE example, on mechanical vibrations, is treated in 
Section 2.2, where we cover many different physical contexts and different 
choices of scales. Scaling the standard, linear wave equation is the topic of 
Chapter 3.1, with discussion of how boundary and initial conditions influence 
the choice of scales. Another PDE example, the diffusion equation, appears 
in Chapter 3.2. Here we progress from a simple linear diffusion equation in 
1D to a study of how scales are influenced by an oscillatory boundary con- 
dition. Nonlinear diffusion models, as well as convection-diffusion PDEs, are 
elaborated on. The final Chapter is devoted to many famous PDEs arising 
from continuum models: elasticity, viscous fluid flow, thermal convection, etc. 

The mathematics is translated into complete computer codes for the ODE 
and simpler PDE problems. 

Experimental fluid mechanics is a field full of relations involving dimen- 
sionless numbers such as the Grashof and Prandtl numbers, but none of the 
textbooks the authors have seen explain how these numbers actually relate to 
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dimensionless forms of the governing equations. Consequently, this non-trivial 
topic is particularly highlighted in the fluid mechanics examples. 

The mathematics in the first two chapters is very gentle and requires 
no more background than basic one-variable calculus and preferably some 
knowledge of differential equation models. The next chapter involves PDEs 
and assumes familiarity with basic models for wave phenomena, diffusion, 
and combined convection-diffusion. The final chapter is meant for readers 
with knowledge of the physics and mathematics of continuum mechanical 
models. The mathematical level of the text rises quickly after the first two 
chapters. 

In the first two chapters, much of the mathematics is accompanied by com- 
plete (yet short) computer codes. The programming level requires familiarity 
with procedural programming in Python. As the mathematical level rises, 
the computer codes get much more comprehensive, and we refer to some files 
for computational examples in chapter three. 

The pedagogy is to saturate the reader with lots of detailed examples to 
provide an understanding for the topic, primarily because the choice of scales 
depends on the problem at hand. One can also view the notes as a reference 
on how to scale many of the most important differential equation models in 
physics. For the simpler differential equations in Chapters 2 and 3, we present 
computer code for many computational examples, but the treatment of the 
advanced models in Chapter 4 is more superficial to limit the size of that 
chapter. 

The exercises are named either Exercise or Problem. The latter is a stand- 
alone exercise without reference to the rest of the text, while the former 
typically extends a topic in the text or refers to sections or formulas in the 
text. 


What this booklet is and is not 


Books containing material on scaling and non-dimensionalization very 
often cover topics not treated in the present notes, e.g., the key topic 
of dimensional analysis and the famous Buckingham Pi Theorem [1, 
8], which we discuss only briefly in section 1.1.3. Similarly, analytical 
solution methods like perturbation techniques and similarity solutions, 
which represent classical methods closely related to scaling and non- 


dimensionalization, are not addressed herein. There are numerous texts 
on perturbation techniques, and these methods build on an already 
scaled differential equations. Similarity solutions do not fit within the 
present scope since these involve non-dimensional combinations of the 
unscaled independent variables to derive new differential equations that 
are easier to solve. 

Our scope is to scale differential equations to simplify the setting of 
parameters in numerical simulations, and at the same time understand 
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more of the physics through interpretation of the dimensionless numbers 
that automatically arise from the scaling procedure. 


With these notes, we hope to demystify the thinking involved in scale 
determination and encourage numerical simulations to be performed with 
dimensionless differential equation models. 

All program and data files referred to in this book are available from the 
book’s primary web site: URL: http://hplgit.github.io/scaling-book/ 
doc/web/. This site also features a version of the book with exercises. 


Acknowledgments. Professor Svein Linge provided very detailed, construc- 
tive comments on the entire manuscript and helped improve the reading qual- 
ity significantly. Yapi Donatien Achou assisted with proof reading. Significant 
portions of the present text were written when the first author was fed with 
FOLFIRINOX (and thereby kept alive) by Linda Falch-Koslung, Dr. Olav 
Dajani, and the rest of the OUS team. There would simply be no booklet 
without their efforts. It is also a great pleasure to express our sincere thanks to 
the Springer and Simula team that handled the prompt editing and produc- 
tion of the text: Martin Peters, Ruth Allewelt, Aslak Tveito, and Asmund 
Ødegård. 


Oslo, November 2015 Hans Petter Langtangen, Geir K. Pedersen 
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Chapter 1 
Dimensions and units 


A mechanical system undergoing one-dimensional damped vibrations can be 
modeled by the equation 


mu” +bu' +ku=0, (1.1) 


where m is the mass of the system, b is some damping coefficient, k is a spring 
constant, and u(t) is the displacement of the system. This is an equation ex- 
pressing the balance of three physical effects: mu” (mass times acceleration), 
bu’ (damping force), and ku (spring force). The different physical quantities, 
such as m, u(t), b, and k, all have different dimensions, measured in different 
units, but mu”, bu’, and ku must all have the same dimension, otherwise it 
would not make sense to add them. 


1.1 Fundamental concepts 
1.1.1 Base units and dimensions 


Base units have the important property that all other units derive from them. 
In the SI system, there are seven such base units and corresponding physical 
quantities: meter (m) for length, kilogram (kg) for mass, second (s) for time, 
kelvin (K) for temperature, ampere (A) for electric current, candela (cd) for 
luminous intensity, and mole (mol) for the amount of substance. 

We need some suitable mathematical notation to calculate with dimensions 
like length, mass, time, and so forth. The dimension of length is written as 
[L], the dimension of mass as [M], the dimension of time as [T], and the 
dimension of temperature as [O] (the dimensions of the other base units 
are simply omitted as we do not make much use of them in this text). The 
dimension of a derived unit like velocity, which is distance (length) divided by 
time, then becomes [LT~'] in this notation. The dimension of force, another 
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derived unit, is the same as the dimension of mass times acceleration, and 
hence the dimension of force is [MLT~?]. 

Let us find the dimensions of the terms in (1.1). A displacement u(t) has 
dimension [L]. The derivative u’(t) is change of displacement, which has di- 
mension [L], divided by a time interval, which has dimension [T], implying 
that the dimension of u’ is [LT~']. This result coincides with the interpre- 
tation of u’ as velocity and the fact that velocity is defined as distance ([L]) 
per time ([T]). 

Looking at (1.1), and interpreting u(t) as displacement, we realize that 
the term mu” (mass times acceleration) has dimension [MLT~?]. The term 
bu’ must have the same dimension, and since u’ has dimension Ss ame b 
must have dimension [MT~"]. Finally, ku must also have dimension [MLT ’?], 
implying that k is a parameter with dimension [MT~?]. 

The unit of a physical quantity follows from the dimension expression. For 
example, since velocity has dimension [LT =t) and length is measured in m 
while time is measured in s, the unit for velocity becomes m/s. Similarly, 
force has dimension [MLT~?] and unit kg m/s”. The k parameter in (1.1) is 


measured in kgs~?. 


Dimension of derivatives 


The easiest way to realize the dimension of a derivative, is to express 
the derivative as a finite difference. For a function u(t) we have 


du u(t + At) — u(t) 

dt ~ At 
where At is a small time interval. If u denotes a velocity, its dimension 
is [LT]~', and u(t+ At) — u(t) gets the same dimension. The time in- 
terval has dimension [T], and consequently, the finite difference gets the 
dimension [LT]~?. In general, the dimension of the derivative du/dt is 
the dimension of u divided by the dimension of t. 


1.1.2 Dimensions of common physical quantities 


Many derived quantities are measured in derived units that have their own 
name. Force is one example: Newton (N) is a derived unit for force, equal 
to kg m/s”. Another derived unit is Pascal (Pa) for pressure and stress, i.e., 
force per area. The unit of Pa then equals N/ m? or kg/ ms”. Below are more 
names for derived quantities, listed with their units. 


1.1 Fundamental concepts 


Name Symbol Physical quantity Unit 
radian rad angle 1 

hertz Hz frequency st 
newton N force, weight kg m/s” 
pascal Pa pressure, stress N/m? 
joule J energy, work, heat Nm 
watt W power J/s 


Some common physical quantities and their dimensions are listed next. 


Quantity 


stress 

pressure 

density 

strain 

Young’s modulus 
Poisson’s ratio 
Lame’ parameters À and yu 
moment (of a force) 
impulse 

linear momentum 
angular momentum 
work 

energy 

power 

heat 

heat flux 
emperature 

heat capacity 
specific heat capacity 
hermal conductivity 
dynamic viscosity 
kinematic viscosity 


surface tension 


Relation 


orce/area 

force/area 

mass/volume 

displacement /length 
stress/strain 

ransverse strain/axial strain 
stress/strain 

distance x force 

orce x time 

mass x velocity 

distance x mass x velocity 
orce x distance 


work 
work/time 
work 

heat rate/area 
base unit 


heat change/temperature change 


heat capacity/unit mass 


heat flux/temperature gradient 
shear stress/velocity gradient 


dynamic viscosity /density 
energy /area 


Unit Dimension 
N/m? = Pa [MT~?L71] 
N/m? = Pa MT~?L7 14] 
kg/m? ML 3} 

1 1] 
N/m? = Pa [MT~?L7}] 
1 1] 
N/m? = Pa [MT~?L7}] 
Nm ML?T~? 
Ns MLT~4] 

kg m/s MLT~ 1} 

kg m?/s  [ML?TT! 
Nm=J ML?T~? 
Nm = J ML?T~? 
Nm/s=W [ML?T~3 

J ML?T~? 
Wm? MT~ 3] 

K ©] 

J/K ML?T-?07}] 
JK~'kg7* [L?T~2071] 
Wm tK! [MLT~367*] 
kgm™ts7} [ML ITT} 
m?/s part] 
J/m? MT~?| 


Prefixes for units. Units often have prefixes!. For example, kilo (k) is a 
prefix for 1000, so kg is 1000 g. Similarly, GPa means giga pascal or 10° Pa. 


1.1.3 The Buckingham Pi theorem 


Almost all texts on scaling has a treatment of the famous Buckingham Pi the- 
orem, which can be used to derive physical laws based on unit compatibility 
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rather than the underlying physical mechanisms. This booklet has its focus 
on models where the physical mechanisms are already expressed through dif- 
ferential equations. Nevertheless, the Pi theorem has a remarkable position 
in the literature on scaling, and since we will occasionally make references to 
it, the theorem is briefly discussed below. 

The theorem itself is simply stated in two parts. First, if a problem in- 
volves n physical parameters in which m independent unit-types (such as 
length, mass etc.) appear, then the parameters can be combined to exactly 
n— m independent dimensionless numbers, referred to as Pi’s. Second, any 
unit-free relation between the original n parameters can be transformed into 
a relation between the n — m dimensionless numbers. Such relations may be 
identities or inequalities stating, for instance, whether or not a given effect is 
negligible. Moreover, the transformation of an equation set into dimension- 
less form corresponds to expressing the coefficients, as well as the free and 
dependent variables, in terms of Pi’s. 

As an example, think of a body moving at constant speed v. What is the 
distance s traveled in time t? The Pi theorem results in one dimensionless 
variable 7 = vt/s and leads to the formula s = Cwt, where C is an undeter- 
mined constant. The result is very close to the well-known formula s = vt 
arising from the differential equation s’ = v in physics, but with an extra 
constant. 

At first glance the Pi theorem may appear as bordering on the trivial. 
However, it may produce remarkable progress for selected problems, such as 
turbulent jets, nuclear blasts, or similarity solutions, without the detailed 
knowledge of mathematical or physical models. Hence, to a novice in scaling 
it may stand out as something very profound, if not magical. Anyhow, as 
one moves on to more complex problems with many parameters, the use of 
the theorem yields comparatively less gain as the number of Pi’s becomes 
large. Many Pi’s may also be recombined in many ways. Thus, good physical 
insight, and/or information conveyed through an equation set, is required 
to pick the useful dimensionless numbers or the appropriate scaling of the 
said equation set. Sometimes scrutiny of the equations also reveals that some 
Pi’s, obtained by applying the theorem, in fact may be removed from the 
problem. As a consequence, when modeling a complex physical problem, the 
real assessment of scaling and dimensionless numbers will anyhow be included 
in the analysis of the governing equations instead of being a separate issue 
left with the Pi theorem. In textbooks and articles alike, the discussion of 
scaling in the context of the equations are too often missing or presented in a 
half-hearted fashion. Hence, the authors’ focus will be on this process, while 
we do not provide much in the way of examples on the Pi theorem. We do 
not allude that the Pi theorem is of little value. In a number of contexts, 
such as in experiments, it may provide valuable and even crucial guidance, 
but in this particular textbook we seek to tell the complementary story on 
scaling. Moreover, as will be shown in this booklet, the dimensionless numbers 
in a problem also arise, in a very natural way, from scaling the differential 
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equations. Provided one has a model based on differential equations, there is 
actually no need for classical dimensional analysis. 


1.1.4 Absolute errors, relative errors, and units 


Mathematically, it does not matter what units we use for a physical quantity. 
However, when we deal with approximations and errors, units are important. 
Suppose we work with a geophysical problem where the length scale is typ- 
ically measured in km and we have an approximation 12.5 km to the exact 
value 12.52 km. The error is then 0.02 km. Switching units to mm leads to 
an error of 20,000 mm. A program working in mm would report 2-10° as the 
error, while a program working in km would print 0.02. The absolute error 
is therefore sensitive to the choice of units. This fact motivates the use of 
relative error: (exact - approximate)/exact, since units then cancel. In the 
present example, one gets a relative error of 1.6-107% regardless of whether 
the length is measured in km or mm. 

Nevertheless, rather than relying solely on relative errors, it is in general 
better to scale the problem such that the quantities entering the computations 
are of unit size (or at least moderate) instead of being very large or very small. 
The techniques of these notes show how this can be done. 


1.1.5 Units and computers 


Traditional numerical computing involves numbers only and therefore re- 
quires dimensionless mathematical expressions. Usually, an implicit trivial 
scaling is used. One can, for example, just scale all length quantities by 1 
m, all time quantities by 1 s, and all mass quantities by 1 kg, to obtain 
the dimensionless numbers needed for calculations. This is the most common 
approach, although it is very seldom explicitly stated. 

Symbolic computing packages, such as Mathematica and Maple, allow 
computations with quantities that have dimension. This is also possible in 
popular computer languages used for numerical computing (Section 1.1.8 
provides a specific example in Python). 


1.1.6 Unit systems 


Confusion arises quickly when some physical quantities are expressed in SI 
units while others are in US or British units. Density could, for instance, be 
given in unit of ounce per teaspoon. Although unit conversion tables are fre- 
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quently met in school, errors in unit conversion probably rank highest among 
all errors committed by scientists and engineers (and when a unit conversion 
error makes an airplane’s fuel run out, it is serious!). Having good soft- 
ware tools to assist in unit conversion is therefore paramount, motivating the 
treatment of this topic in Sections 1.1.8 and 1.2. Readers who are primarily 
interested in the mathematical scaling technique may safely skip this material 
and jump right to Section 2.1. 


1.1.7 Example on challenges arising from unit systems 


A slightly elaborated example on scaling in an actual science/engineering 
project may stimulate the reader’s motivation. In its full extent, the study 
of tsunamis spans geophysics, geology, history, fluid dynamics, statistics, 
geodesy, engineering, and civil protection. This complexity reflects in a diver- 
sity of practices concerning the use of units, scales, and concepts. If we narrow 
the scope to modeling of tsunami propagation, the scaling aspect, at least, 
may seem simple as we are mainly concerned with length and time. Still, even 
here the non-uniformity concerning physical units is an encumbrance. 

A minor issue is the occasional use of non-SI units such as inches, or in old 
charts, even fathoms. More important is the non-uniformity in the magnitude 
of the different variables, and the differences in the inherent horizontal and 
vertical scales in particular. Typically, surface elevations are in meters or 
smaller. For far-field deep water propagation, as well as small tsunamis (which 
are still of scientific interest) surface elevations are often given in cm or even 
mm. In the deep ocean, the characteristic depth is orders of magnitude larger 
than this, typically 5000m. Propagation distances, on the other hand, are 
hundreds or thousands of kilometers. Often locations and computational grids 
are best described in geographical coordinates (longitude/latitude) which are 
related to SI units by 1 latitude minute being roughly one nautical mile 
(1852m), and 1 longitude minute being this quantity times the cosine of the 
latitude. Wave periods of tsunamis mostly range from minutes to an hour, 
hopefully sufficiently short to be well separated from the half-daily period of 
the tides. Propagation times are typically hours or maybe the better part of 
a day when the Pacific Ocean is traversed. 

The scientists, engineers, and bureaucrats in the tsunami community tend 
to be particular and non-conform concerning formats and units, as well as 
the type of data required. To accommodate these demands, a tsunami mod- 
eler must produce a diversity of data which are in units and formats which 
cannot be used internally in her models. On the other hand, she must also be 
prepared to accept the input data in diversified forms. Some data sets may 
be large, implying that unnecessary duplication, with different units or scal- 
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ing, should be avoided. In addition, tsunami models are often bench-marked 
through comparison with experimental data. The lab scale is generally cm 
or m, at most, which implies that measured data are provided in different 
units (than used in real earth-scale events), or even in volts, with conversion 
information, as obtained from the measuring gauges. 

All the unit particulars in various file formats is clearly a nuisance and give 
rise to a number of misconceptions and errors that may cause loss of precious 
time or efforts. To reduce such problems, developers of computational tools 
should combine a reasonable flexibility concerning units in input and output 
with a clear and consistent convention for scaling within the tools. In fact, 
this also applies to academic tools for in-house use. 

The discussion above points to some best practices that these notes pro- 
motes. First, always compute with scaled differential equation models. This 
booklet tells you how to do that. Second, users of software often want to 
specify input data with dimension and get output data with dimension. The 
software should then apply tools like PhysicalQuantity (Section 1.1.8) or 
the more sophisticated Parampool package (Section 1.2) to allow input with 
explicit dimensions and convert the dimensions to the right types if necessary. 
It is trivial to apply these tools if the computational software is written in 
Python, but it is even straightforward if the software is written in compiled 
languages like Fortran, C, or C++. In the latter case one just makes an input 
reading module in Python that grabs data from a user interface and feeds 
them into the computational software, either through files or function calls 
(the relevant functions to be called must be wrapped in Python with tools 
like f2py?, Cython*, Weave®, SWIG®, Instant”, or similar, see [7, Appendix 
C] for basic examples on f2py and Cython wrapping of C and Fortran code). 


1.1.8 PhysicalQuantity: a tool for computing with units 


These notes contain quite some computer code to illustrate how the theory 
maps in detail to running software. Python is the programming language 
used, primarily because it is an easy-to-read, powerful, full-fledged language 
that allows MATLAB-like code as well as class-based code typically used in 
Java, C#, and C++. The Python ecosystem for scientific computing has in 
recent years grown fast in popularity and acts as a replacement for more spe- 
cialized tools like MATLAB, R, and IDL. The coding examples in this booklet 
requires only familiarity with basic procedural programming in Python. 


3http://docs.scipy .org/doc/numpy-dev/f2py/ 

4nttp://cython.org/ 
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Readers without knowledge of Python variables, functions, if tests, and 
module import should consult, e.g., a brief tutorial on scientific Python®, 
the Python Scientific Lecture Notes, or a full textbook [4] in parallel with 
reading about Python code in the present notes. 


These notes apply Python 2.7 


Python exists in two incompatible versions, numbered 2 and 3. The 
differences can be made small, and there are tools to write code that 
runs under both versions. 


As Python version 2 is still dominating in scientific computing, we 
stick to this version, but write code in version 2.7 that is as close as 
possible to version 3.4 and later. In most of our programs, only the 
print statement differs between version 2 and 3. 


Computations with units in Python are well supported by the very use- 
ful tool PhysicalQuantity from the ScientificPython package!? by Kon- 
rad Hinsen. Unfortunately, ScientificPython does not, at the time of this 
writing, work with NumPy version 1.9 or later, so we have isolated the 
PhysicalQuantity object in a module PhysicalQuantities!! and made 
it publicly available on GitHub. There is also an alternative package Unum! 
for computing with numbers with units, but we shall stick to the former 
module here. 

Let us demonstrate the usage of the PhysicalQuantity object by com- 
puting s = vt, where v is a velocity given in the unit yards per minute and t 
is time measured in hours. First we need to know what the units are called 
in PhysicalQuantities. To this end, run pydoc PhysicalQuantities, or 


Terminal 


Terminal> pydoc Scientific.Physics.PhysicalQuantities 


if you have the entire ScientificPython package installed. The resulting docu- 
mentation shows the names of the units. In particular, yards are specified by 
yd, minutes by min, and hours by h. We can now compute s = vt as follows: 


>>> # With ScientificPython: 
>>> from Scientific.Physics.PhysicalQuantities import \ 

. PhysicalQuantity as PQ 
>>> # With PhysicalQuantities as separate/stand-alone module: 
>>> from PhysicalQuantities import PhysicalQuantity as PQ 
>>> 


Shttp://hplgit.github.io/bumpy/doc/web/index .html 
°nttp://scipy-lectures.github.com/ 
10nttps://bitbucket .org/khinsen/scientificpython 
llhttps://github.com/hplgit/physical-quantities 
2nttps://bitbucket .org/kiv/unum/ 
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>>> v = PQ(’?120 yd/min’) # velocity 


>>> t = PQ(’1 h’) # time 
>>> s = vet # distance 
>>> print s # s is string 


120.0 h*yd/min 


The odd unit h*yd/min is better converted to a standard SI unit such as 
meter: 


>>> s.convertToUnit(’m’) 
>>> print s 
6583.68 m 


Note that s is a PhysicalQuantity object with a value and a unit. For 
mathematical computations we need to extract the value as a float object. 
We can also extract the unit as a string: 


>>> print s.getValue() # float 
6583.68 

>>> print s.getUnitName() # string 
m 


Here is an example on how to convert the odd velocity unit yards per 
minute to something more standard: 


>>> v.convertToUnit(’km/h’) 
>>> print v 

6.58368 km/h 

>>> v.convertToUnit(’m/s’) 
>>> print v 

1.8288 m/s 


As another example on unit conversion, say you look up the specific heat 
capacity of water to be 1 calg~!K~!. What is the corresponding value in the 
standard unit Jg~'K~' where joule replaces calorie? 


>>> c = PQ(’1 cal/(g*K)’) 

>>> c.convertToUnit (’ J/(g*K)’) 
>>> print c 

4.184 J/K/g 


1.2 Parampool: user interfaces with automatic 
unit conversion 


The Parampool!? package allows creation of user interfaces with support 
for units and unit conversion. Values of parameters can be set as a number 
with a unit. The parameters can be registered beforehand with a preferred 


13https://github.com/hplgit/parampool 
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unit, and whatever the user prescribes, the value and unit are converted so 
the unit becomes the registered unit. Parampool supports various type of 
user interfaces: command-line arguments (option-value pairs), text files, and 
interactive web pages. All of these are described next. 


Example application. As case, we want to make software for computing 
with the simple formula s = vot + Zat. We want vo to be a velocity with unit 


m/s, a to be acceleration with unit m/s?, t to be time measured in s, and 
consequently s will be a distance measured in m. 


1.2.1 Pool of parameters 


First, Parampool requires us to define a pool of all input parameters, which 
is here simply represented by list of dictionaries, where each dictionary holds 
information about one parameter. It is possible to organize input parameters 
in a tree structure with subpools that themselves may have subpools, but 
for our simple application we just need a flat structure with three input 
parameters: vo, a, and t. These parameters are put in a subpool called “Main”. 
The pool is created by the code 


def define_input(): 
pool = [ 

‘Maan? [ 
dict (name=’initial velocity’, default=1.0, unit=’m/s’), 
dict (name=’acceleration’, default=1.0, unit=’m/s**2’), 
dict (name=’time’, default=10.0, unit=’s’) 
] 

] 


from parampool.pool.UI import listtree2Pool 
pool = listtree2Pool(pool) # convert list to Pool object 
return pool 


For each parameter we can define a logical name, such as initial velocity, 
a default value, and a unit. Additional properties are also allowed, see the 


Parampool documentation’. 


Tip: specify default values of numbers as float objects 


Note that we do not just write 1, but 1.0 as default. Had 1 been used, 
Parampool would have interpreted our parameter as an integer and 


would therefore convert input like 2.5 m/s to 2 m/s. To ensure that a 
real-valued parameter becomes a float object inside the pool, we must 
specify the default value as a real number: 1. or 1.0. (The type of 


M4pttp://hplgit.github.io/parampool/doc/web/index. html 
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an input parameter can alternatively be set explicitly by the str2type 
property, e.g., str2type=float.) 


1.2.2 Fetching pool data for computing 


We can make a little function for fetching values from the pool and computing 
8: 
def distance (pool): 
v_0 = pool.get_value(’initial velocity’) 
a = pool.get_value(’acceleration’) 
t = pool.get_value(’time’) 
s = v_O*t + 0.5%*a*t**x2 
return s 


The pool.get_value function returns the numerical value of the named pa- 
rameter, after the unit has been converted from what the user has specified 
to what was registered in the pool. For example, if the user provides the 
command-line argument -time ’2 h’, Parampool will convert this quantity 
to seconds and pool.get_value(’time’) will return 7200. 


1.2.3 Reading command-line options 


To run the computations, we define the pool, load values from the command 
line, and call distance: 


pool = define_input () 
from parampool.menu.UI import set_values_from_command_line 
pool = set_values_from_command_line (pool) 


s = distance(pool) 
print ’s=/g’ % s 
Parameter names with whitespace must use an underscore for whitespace 
in the command-line option, such as in --Initial_velocity. We can now 
run 


Terminal 


Terminal> python distance.py --initial_velocity ’10 km/h’ \ 
--acceleration 0 --time ’1h 
s=10000 


Notice from the answer (s) that 10 km/h gets converted to m/s and 1 h to s. 
It is also possible to fetch parameter values as PhysicalQuantity objects 
from the pool by calling 
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v_0 = pool.get_value_unit(’Initial velocity’) 


The following variant of the distance function computes with values and 
units: 


def distance_unit(pool): 
# Compute with units 
from parampool.PhysicalQuantities import PhysicalQuantity as PQ 
v_0 = pool.get_value_unit(’initial velocity’) 
a = pool.get_value_unit(’acceleration’) 
t = pool.get_value_unit(’time’) 
s = v_O*t + 0.5*a*t**2 
return s.getValue(), s.getUnitName() 


We can then do 


s, s_unit = distance_unit (pool) 
print ’s=/g’ % s, s_unit 


and get output with the right unit as well. 


1.2.4 Setting default values in a file 


In large applications with lots of input parameters one will often like to define 
a (huge) set of default values specific for a case and then override a few of 
them on the command-line. Such sets of default values can be set in a file 
using syntax like 

subpool Main 

initial velocity = 100 ! yd/min 

acceleration = 0 ! m/s**2 # drop acceleration 

end 


The unit can be given after the ! symbol (and before the comment symbol 
#). 
To read such files we have to add the lines 


from parampool.pool.UI import set_defaults_from_file 
pool = set_defaults_from_file(pool) 


before the call to set_defaults_from_command_line. 

If the above commands are stored in a file distance.dat, we give this file 
information to the program through the option -poolfile distance.dat. 
Running just 


Terminal 


Terminal> python distance.py --poolfile distance.dat 
s=15.25 m 
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first loads the velocity 100 yd/min converted to 1.524 m/s and zero accel- 
eration into the pool system and, and then we call distance_unit, which 
loads these values from the pool along with the default value for time, set as 
10 s. The calculation is then s = 1.524-10+0 = 15.24 with unit m. We can 
override the time and/or the other two parameters on the command line: 


Terminal 


Terminal> python distance.py --poolfile distance.dat --time ’2 h’ 
s=10972.8 m 


The resulting calculations are s = 1.524- 7200 +0 = 10972.8. You are encour- 
aged to play around with the distance.py program. 


1.2.5 Specifying multiple values of input parameters 


Parampool has an interesting feature: multiple values can be assigned to an 
input parameter, thereby making it easy for an application to run through all 
combinations of all parameters. We can demonstrate this feature by making 
a table of vg, a, t, and s values. In the compute function, we need to call 
pool.get_values instead of pool. get_value to get a list of all the values 
that were specified for the parameter in question. By nesting loops over all 
parameters, we visit all combinations of all parameters as specified by the 
user: 
def distance_table(pool): 

"""Grab multiple values of parameters from the pool.""" 

table = [] 

for v_0 in pool.get_values(’initial velocity’): 

for a in pool.get_values(’acceleration’): 
for t in pool.get_values(’time’): 
s = v_O*#t + 0.5*xa*t**2 
table.append((v_0, a, t, s)) 

return table 
In case just a single value was specified for a parameter, pool.get_values 
returns this value only and there will be only one pass in the associated loop. 

After loading command-line arguments into our pool object, we can call 

distance_table instead of distance or distance_unit and write out a 
nicely formatted table of results: 


table = distance_table(pool) 


fom vos a.) e S intable: 
eaaa 2 aek || Os || o || aE E 7 GeO, a, Eo E 
Prin EN SSeS sSS Sse s SSS SSS SSS Se SS SS SSS SS SSS lk 


Here is a sample run, 


14 1 Dimensions and units 


Terminal 


Terminal> python distance.py --time °1h&2h& 3h’ \ 
--acceleration ’0 m/s**2 & 1 m/s**2 & 1 yd/s**2’ \ 
--initial_velocity ’1 & 5’ 


| v_0 | a | t | s | 
| pers rt ss st | 
| 1.000 | 0.000 | 3600.000 | 3600.000 | 
| 1.000 | 0.000 | 7200.000 | 7200.000 | 
| 1.000 | 0.000 | 10800.000 | 10800.000 | 
| 1.000 | 1.000 | 3600.000 | 6483600.000 | 
| 1.000 | 1.000 | 7200.000 | 25927200.000 | 
| 1.000 | 1.000 | 10800.000 | 58330800.000 | 
| 1.000 | 0.914 | 3600.000 | 5928912.000 | 
| 1.000 | 0.914 | 7200.000 | 23708448.000 | 
| 1.000 | 0.914 | 10800.000 | 53338608.000 | 
| 5.000 | 0.000 | 3600.000 | 18000.000 | 
| 5.000 | 0.000 | 7200.000 | 36000.000 | 
| 5.000 | 0.000 | 10800.000 | 54000.000 | 
| 5.000 | 1.000 | 3600.000 | 6498000.000 | 
| 5.000 | 1.000 | 7200.000 | 25956000.000 | 
| 5.000 | 1.000 | 10800.000 | 58374000.000 | 
| 5.000 | 0.914 | 3600.000 | 5943312.000 | 
| 5.000 | 0.914 | 7200.000 | 23737248.000 | 
| 5.000 | 0.914 | 10800.000 | 53381808.000 | 


Notice that some of the multiple values have dimensions different from the 
registered dimension for that parameter, and the table shows that conversion 
to the right dimension has taken place. 


1.2.6 Generating a graphical user interface 


For the fun of it, we can easily generate a graphical user interface via Param- 
pool. We wrap the distance_unit function in a function that returns the 
result in some nice-looking HTML code: 


def distance_unit2(pool): 
# Wrap result from distance_unit in HTML 
s, s_unit = distance_unit (pool) 
return? <p>Distancen </b>. 2i Ase (S Seunat) 


In addition, we must make a file generate_distance_GUI.py with the simple 
content 


from parampool.generator.flask import generate 
from distance import distance_unit2, define_input 


generate(distance_unit2, pool_function=define_input, MathJax=True) 
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Running generate_distance_GUI.py creates a Flask-based web interface! 
to our distance_unit function, see Figure 1.1. The text fields in this GUI 
allow specification of parameters with numbers and units, e.g., acceleration 
with unit yards per minute squared, as shown in the figure. Hovering the 
mouse slightly to the left of the text field causes a little black window to pop 
up with the registered unit of that parameter. 


Input: Result: 
open all | close all Distance: 200592.00 m 
Main 
initial velocity 10 
acceleration 1.0 yd/min**2 
time 10.0h 
Compute 


Fig. 1.1 Web GUI where parameters can be specified with units. 


With examples shown above, the reader should be able to make use of 
the PhysicalQuantity object and the Parampool package in programs and 
thereby work safely with units. For the coming text, where we discuss the 
craft of scaling in detail, we shall just work in standard SI units and avoid unit 
conversion so there will be no more use of PhysicalQuantity and Parampool. 


Open Access This chapter is distributed under the terms of the Creative Commons Attribution- 
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15You need to have Flask and additional packages installed. This is easy to do with 
a few pip install commands, see [5] or [6]. 


Chapter 2 
Ordinary differential equation models 


This chapter introduces the basic techniques of scaling and the ways to reason 
about scales. The first class of examples targets exponential decay models, 
starting with the simple ordinary differential equation (ODE) for exponential 
decay processes: u’ = —au, with constant a > 0. Then we progress to vari- 
ous generalizations of this ODE, including nonlinear versions and systems of 
ODEs. The next class of examples concerns second-order ODEs for oscilla- 
tory systems, where the simplest ODE reads mu” + ku = 0, with m and k 
as positive constants. Various extensions with damping and force terms are 
discussed in detail. 


2.1 Exponential decay problems 


2.1.1 Fundamental ideas of scaling 


Scaling is an extremely useful technique in mathematical modeling and nu- 
merical simulation. The purpose of the technique is three-fold: 


1. Make independent and dependent variables dimensionless. 
2. Make the size of independent and dependent variables about unity. 
3. Reduce the number of independent physical parameters in the model. 


The first two items mean that for any variable, denote it by q, we introduce 
a corresponding dimensionless variable 
oo GH go 
q= , 
qe 
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where qo is a reference value of q (qo = 0 is a common choice) and qe is a 
characteristic size of |q|, often referred to as “a scale”. Since the numerator and 
denominator have the same dimension, q becomes a dimensionless number. 

If qc is the maximum value of |q — qo|, we see that 0 < |g| < 1. How to find 
qe is sometimes the big challenge of scaling. Examples will illustrate various 
approaches to meet this challenge. 

The many coming examples on scaling differential equations contain the 
following pedagogical ingredients to meet the desired learning outcomes. 


e Teach the technical steps of making a mathematical model, based on dif- 
ferential equations, dimensionless. 

e Describe various techniques for reasoning about the scales, i.e., finding the 
characteristic sizes of quantities. 

e Teach how to identify and interpret dimensionless numbers arising from 
the scaling process. 

e Provide a lot of different examples on making models dimensionless with 
physically correct scales. 

e Show how symbolic software (SymPy) can be used to derive exact solutions 
of differential equations. 

e Explain how to run a dimensionless model with software developed for the 
problem with dimensions. 


2.1.2 The basic model problem 


Processes undergoing exponential reduction can be modeled by the ODE 
problem 


u(t) =—au(t), u(0) =T, (2.1) 


where a, I > 0 are prescribed parameters, and u(t) is the unknown function. 
For the particular model with a constant a, we can easily derive the exact 
solution, u(t) = Ie~**, which is helpful to have in mind during the scaling 
process. 


Example: Population dynamics. The evolution of a population of hu- 
mans, animals, cells, etc., under unlimited access to resources, can be mod- 
eled by (2.1). Then u is the number of individuals in the population, strictly 
speaking an integer, but well modeled by a real number in large populations. 
The parameter a is the increase in the number of individuals per time and 
per individual. 


Example: Decay of pressure with altitude. The simple model (2.1) also 
governs the pressure in the atmosphere (under many assumptions, such air is 
an ideal gas in equilibrium). In this case u is the pressure, measured in Nm7?; 
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t is the height in meters; and a= _M/(R*T), where M is the molar mass of the 
Earth’s air (0.029 kg/mol), R* is the universal gas constant (8.314 At), 
and T is the temperature in Kelvin (K). The temperature depends on the 


height so we have a = a(t). 


2.1.3 The technical steps of the scaling procedure 


Step 1: Identify independent and dependent variables. There is one 
independent variable, t, and one dependent variable, u. 


Step 2: Make independent and dependent variables dimensionless. 
We introduce a new dimensionless t, called t, defined by 
a et 
t= —, (2.2) 
te 
where te is a characteristic value of t. Similarly, we introduce a dimensionless 
u, named u, according to 
“a=—, (2.3) 
Uc 
where uc is a constant characteristic size of u. When wu has a specific inter- 
pretation, say when (2.1) models pressure in an atmospheric layer, ue would 
be referred to as characteristic pressure. For a decaying population, ue may 
be a characteristic number of members in the population, e.g., the initial 
population I. 


Step 3: Derive the model involving only dimensionless variables. 
The next task is to insert the new dimensionless variables in the governing 
mathematical model. That is, we replace t by tet and u by ucù in (2.1). The 
derivative with respect to t is derived through the chain rule as 


du d(ucu) dt dul ucdu 


dt dt d "dite te dt’ 
The model (2.1) now becomes 
Uc du 7 P 
oa =—aucū, ucU0)=T. (2.4) 


Step 4: Make each term dimensionless. Equation (2.4) still has terms 
with dimensions. To make each term dimensionless, we usually divide by the 
coefficient in front of the term with the highest time derivative (but dividing 
by any coefficient in any term will do). The result is 


du 


“=at, (0) = uz". (2.5) 
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Step 5: Estimate the scales. A characteristic quantity like te reflects the 
time scale in the problem. Estimating such a time scale is certainly the most 
challenging part of the scaling procedure. There are different ways to reason. 
The first approach is to aim at a size of u and its derivatives that is of order 
unity. If ue is chosen such that |u| is of size unity, we see from (2.5) that 
du/dt is of the size of U (i.e., unity) if we choose te = 1/a. 

Alternatively, we may look at a special case of the model where we have 
analytical insight that can guide the choice of scales. In the present problem 
we are lucky to know the exact solution for any value of the input data as 
long as a is a constant. For exponential decay, u(t) ~ e~%, it is common to 
define a characteristic time scale te as the time it takes to reduce the initial 
value of u by a factor of 1/e (also called the e-folding time): 


= Iri = = 
e tte — Le a-0 = e ate — 1 
e 
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from which it follows that te =1/a. Note that using an exact solution of the 
problem to determine scales is not a requirement, just a useful help in the 
few cases where we actually have access to an exact solution. 

In this example, two different, yet common ways of reasoning, lead to the 
same value of te. However, instead of using the e-folding time we could use 
the half-time of the exponential decay as characteristic time, which is also a 
very common measure of the time scale in such processes. The half time is 
defined as the time it takes to halve u: 


1 
Uae ere => te=a'In2. 


There is a factor ln 2 = 0.69 difference from the other te value. As long as the 
factor is not an order of magnitude or more different, we do not pay attention 
factors like In2 and skip them, simply to make formulas look nicer. Using 


te =a~'In2 as time scale leads to a scaled differential equation u’ = —(In2)u, 
which is fine, but an unusual form. People tend to prefer the simpler ODE 
u’ = —u, which arises from te = 1/a, and we shall therefore use this time scale. 


Regarding uc, we may look at the initial condition and realize that the 
choice ue = I makes u(0) = 1. For t > 0, the differential equation expresses 
explicitly that u decreases, so uc = I gives u € (0, 1]. Scaling a variable q such 
that |q| € [0,1] is always the ultimate goal, and this goal is in fact obtained 
here! Next best result is to ensure that the magnitude of |q| is not “big” or 
“small”, in the sense that the size is neither as large as 10 or 100, nor as small 
as 0.1 or 0.01. (In the present problem, where we are lucky to have an exact 
solution u(t) = Ie~%, we may look at this to explicitly see that u € (0,/] such 
that ue = I gives ù € (0,1]). 

With te =1/a and uc = I, we have the final dimensionless model 


dū 


=, a(0)=1. (2.6) 
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This is a remarkable result in the sense that all physical parameters (a and I) 
are removed from the model! Or more precisely, there are no physical input 
parameters to assign before using the model. In particular, numerical inves- 
tigations of the original model (2.1) would need experiments with different a 
and I values, while numerical investigations of (2.6) can be limited to a single 


run! As soon as we have computed the curve u(t), we can find the solution 
u(t) of (2.1) by 


u(t) = ucii(t/te) = Tii(at). (2.7) 


This particular transformation actually means stretching the t and wu axes in 
a plot of u(t) by the factors a and I, respectively. 
It is very common to drop the bars when the scaled problem has been 


derived and work further with (2.6) simply written as 


d 
a =-u, u(0)=1. 
Nevertheless, in this booklet we have decided to stick to bars for all dimen- 
sionless quantities. 


2.1.4 Making software for utilizing the scaled model 


Software for solving (2.1) could take advantage of the fact that only one 
simulation of (2.6) is necessary. As soon as we have w(t) accessible, a simple 
scaling (2.7) computes the real u(t) for any given input data a and J. Although 
the numerical computation of u(t) from (2.1) is very fast in this simple model 
problem, using (2.7) is very much faster. In general, a simple rescaling of a 
scaled solution is extremely more computationally efficient than solving a 
differential equation problem. 

We can compute with the dimensionless model (2.6) in two ways, either 
make a solver for (2.6), or reuse a solver for (2.1) with J = 1 and a = 1. We will 
choose the latter approach since it has the advantage of giving us software 
that works both with a dimensionless model and a model with dimensions 
(and all the original physical parameters). 


Software for the original unscaled problem. Assume that we have some 
module decay.py that offers the following functions: 


e solver(I, a, T, dt, theta=0.5) for returning the solution arrays u 
and t, over a time interval [0,7], for (2.1) solved by the so-called 0 rule. 
This rule includes the Forward Euler scheme (0 = 0), the Backward Euler 
scheme (6 = 1), or the Crank-Nicolson (centered midpoint) scheme (6 = $). 

e read_command_line_argparse() for reading parameters in the problem 
from the command line and returning them: I, a, T, theta (0), and a list 


of At values for time steps. (We shall only make use of the first At value.) 
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The basic statements for solving (2.1) are then 


from decay import solver, read_command_line_argparse 
I, a, T, theta, dt_values = read_command_line_argparse() 
u, t = solver(I, a, T, dt_values[0], theta) 


from matplotlib.pyplot import plot, show 
plot(t, u) 
show() 


The module decay .py is developed and explained in Section 5.1.7 in [3]. 

_ To solve the dimensionless problem, just fix I = 1 and a = 1, and choose 
T and At: 

_, _» T, theta, dt_values = read_command_line_argparse() 

u, t = solver(I=1, a=1, T=T, dt=dt_values[0], theta=theta) 


The first two variables returned from read_command_line_argparse are I 
and a, which are ignored here. To indicate that these variables are not to be 
used, we use a “dummy name”, often taken to be the underscore symbol in 
Python. The user can set -I and -a on the command line, since the decay 
module allows this, but we hope the code above has a form that reminds the 
user that these options are not to be used. Also note that T and dt_values [0] 
set on the command line are the desired parameters for solving the scaled 
problem. 


Software for the scaled problem. Turning now to the scaled problem, 
the solver function (originally designed for the unscaled problem) will be 
reused, but it will only be run if it is strictly necessary. That is, when the 
user requests a solution, our code should first check whether that solution can 
be provided by simply scaling a solution already computed and available in a 
file. If not, we will compute an appropriate scaled solution, find the requested 
unscaled solution for the user, and also save the new scaled solution to file 
for possible later use. 

A very plain solution to the problem is found in the file decay_scaled_ 
vi.py. The np.savetxt function saves a two-dimensional array (“table”) to 
a text file, and the np.loadtxt function can load the data back into the 
program. A better solution to this problem is obtained by using the joblib 
package as described next. 


Implementation with joblib. The Python package joblib has function- 
ality that is very convenient for implementing the solver_scaled function. 
The first time a function is called with a set of arguments, the statements in 
the function are executed and the return value is saved to file. If the function 
is called again with the same set of arguments, the statements in the func- 
tion are not executed, but the return value is read from file (of course, many 
files may be stored, one for each combination of parameter values). In com- 
puter science, one would say that joblib in this way provides memorization 
functionality for Python functions. This functionality is particularly aimed at 
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large-scale computations with arrays that one would hesitate to recompute. 

We illustrate the technique here in a very simple mathematical context. 
First we make a solver_scaled function for the scaled model that just 

calls up a solver_unscaled (with J =a=1) for the problem with dimensions: 


from decay import solver as solver_unscaled 
import numpy as np 
import matplotlib.pyplot as plt 


def solver_scaled(T, dt, theta): 


Solve u’=-u, u(0)=1 for (0,T] with step dt and theta method. 
print ’Computing the numerical solution’ 
return solver_unscaled(I=1, a=1, T=T, dt=dt, theta=theta) 


Then we create some “computer memory on disk”, i.e., some disk space to 
store the result of a call to the solver_scaled function. Thereafter, we rede- 
fine the name solver_scaled to a new function, created by joblib, which 
calls our original solver_scaled function if necessary and otherwise loads 
data from file: 


import joblib 
disk_memory = joblib.Memory(cachedir=’ temp’) 
solver_scaled = disk_memory.cache(solver_scaled) 


The solutions are actually stored in files in the cache directory temp. 

A typical use case is to read values from the command line, solve the scaled 
problem (if necessary), unscale the solution, and visualize the solution with 
dimension: 


def unscale(u_scaled, t_scaled, I, a): 
return I*u_scaled, a*t_scaled 


from decay import read_command_line_argparse 


def main(): 
# Read unscaled parameters, solve and plot 
I, a, T, theta, dt_values = read_command_line_argparse() 
dt = dt_values[0] # use only the first dt value 
T_bar = a*T 
dt_bar = a*dt 
u_scaled, t_scaled = solver_scaled(T_bar, dt_bar, theta) 
u, t = unscale(u_scaled, t_scaled, I, a) 


plt.figure() 

plt.plot(t_scaled, u_scaled) 

plt.xlabel(’scaled time’); plt.ylabel(’scaled velocity’) 
plt.title(’Universial solution of scaled problem’) 
plt.savefig(’tmpi.png’); plt.savefig(’tmp1.pdf’) 


plt.figure() 

plt.plot(t, u) 

plt.xlabel(’t’); plt.ylabel(’u’) 

plt.title(’I=%g, a=/g, theta=/g’ % (I, a, theta)) 
plt.savefig(’tmp2.png’); plt.savefig(’tmp2.pdf’) 
plt.show() 
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The complete code resides in the file decay_scaled.py. Note from the code 
above that read_command_line_argparse is supposed to read parameters 
with dimensions (but technically, we solve the scaled problem, if strictly nec- 
essary, and unscale the solution). Let us run 


Terminal 


Terminal> python decay_scaled.py --I 8 --a 0.1 --dt 0.01 --T 50 


A plot of the scaled and unscaled solution appears in Figure 2.1. 


Universial solution of scaled problem l=8, a=0.1, theta=0.5 


scaled velocity 
u 


0 1 2 3 4 5 0.0 0.1 0.2 0.3 0.4 0.5 
scaled time t 


Fig. 2.1 Scaled (left) and unscaled (right) exponential decay. 


Note that we write a message Computing the numerical solution in- 
side the solver_scaled function. We can then easily detect when the solu- 
tion is actually computed from scratch and when it is simply read from file 
(followed by the unscaling procedure). Here is a demo: 


Terminal 


Terminal> # Very first run 

Terminal> python decay_scaled.py --T 7 --a 1 --I 0.5 --dt 0.2 
[Memory] Calling __main__--home-hpl... 
solver_scaled-alias(7.0, 0.2, 0.5) 

Computing the numerical solution 


Terminal> # No change of T, dt, theta - can reuse solution in file 
Terminal> python decay_scaled.py --T 7 --a 4 --I 2.5 --dt 0.2 


Terminal> # Change of dt, must recompute 

Terminal> python decay_scaled.py --T 7 --a 4 --I 2.0 --dt 0.5 
[Memory] Calling __main__--home-hpl... 
solver_scaled-alias(7.0, 0.5, 0.5) 

Computing the numerical solution 


Terminal> # Change of dt again, but dt=0.2 is already in a file 
Terminal> python decay_scaled.py --T 7 --a 0.5 --I 1 --dt 0.2 


We realize that joblib has access to all previous runs and does not re- 
compute unless it is strictly required. Our previous implementation without 
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joblib (in decay_scaled_v1i.py) used only one file (for one numerical case) 
and will therefore perform many more calls to solver_unscaled. 


On the implementation of a simple memoize function 


A memoized function recalls previous results when the same set of argu- 
ments is encountered. That is, the function caches its results. A simple 
implementation stores the arguments in a function call and the returned 
results in a dictionary, and if the arguments are seen again, one looks 
up in the dictionary and returns previously computed results: 
class Memoize: 
det Minit Eeli T DE 
self.f = f 
self.memo = {} # map arguments to results 


def __call__(self, *args): 
if not args in self.memo: 
self .memo [args] = self.f(*args) 
return self .memo [args] 


# Wrap my_compute_function(arg1, arg2, ...) 
my_compute_function = Memoize(my_compute_function) 


The memoize functionality in joblib.Memory is more sophisticated and 
can work very efficiently with large array data structures as arguments. 
Note that the simple version above can only be used when all arguments 
to the function f are immutable (since the key in a dictionary has to 
be immutable). 


2.1.5 Scaling a generalized problem 


Now we consider an extension of the exponential decay ODE to the form 


u'(t) = —au(t)+b, u(0)=T. (2.8) 


One particular model, with constant a and b, is a spherical small-sized or- 


ganism falling in air, 
3rdu ( 0 
u = u+ 1), 2.9 
aV I No eel 


where d, H, 0p, 0, V, and g are physical parameters. The function u(t) rep- 
resents the vertical velocity, being positive upwards. We shall use this model 
in the following. 
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Exact solution. It can be handy to have the exact solution for reference, 
in case of constant a and b: 


Solving differential equations in SymPy 


It can be very useful to use a symbolic computation tool such as SymPy 
to aid us in solving differential equations. Let us therefore demonstrate 
how SymPy can be used to find this solution. First we define the pa- 
rameters in the problem as symbols and u(t) as a function: 


>>> from sympy import * 
>>> t, a, b, I = symbols(’t a b I’, real=True, positive=True) 
>>> u = symbols(’u’, cls=Function) 


The next task is to define the differential equation, either as a symbolic 
expression that is to equal zero, or as an equation Eq(lhs, rhs) with 
lhs and rhs as expressions for the left- and right-hand side): 


>>> # Define differential equation 

>>> eq = diff(u(t), t) + atu(t) - b 
>>> # or 

>>> eq = Eq(diff(u(t), t), -a*u(t) + b) 


The differential equation can be solved by the dsolve function, yielding 
an equation of the form u(t) == expression. We want to grab the 
expression on the right-hand side as our solution: 


>>> sol = dsolve(eq, u(t)) 

>>> print sol 

u(t) == (b + exp(a*(C1 - t)))/a 

>>> u = sol.rhs # grab solution 
>>> print u 

(b + exp(a*(C1 - t)))/a 


The solution contains the unknown integration constant C1, which must 
be determined by the initial condition. We form the equation arising 
from the initial condition u(0) = I: 


>>> C1 = symbols(’C1’) 

>>> eq = Eq(u.subs(t, 0), I) # substitute t by 0 in u 
>>> sol = solve(eq, C1) 

>>> print sol 

[log(I*a - b)/a] 


The one solution that was found (stored in a list!) must then be sub- 
stituted back in the expression u to yield the final solution: 
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>>> u = u.subs(C1, sol[0]) 
>>> print u 
(b + exp(a*(-t + log(I*a - b)/a)))/a 


As in mathematics with pen and paper, we strive to simplify expressions 
also in symbolic computing software. This frequently requires some trial 
and error process with SymPy’s simplification functions. A very stan- 


dard first try is to expand everything and run simplification algorithms: 


>>> u = simplify (expand(u) ) 
>>> print u 
(Ixa + b*exp(a*t) - b)*exp(-ax*t)/a 


Doing latex(u) automatically converts the expression to ATEX syntax 
for inclusion in reports. 


The reader may wonder why we bother with scaling of differential equa- 
tions if SymPy can solved the problem in a nice, closed formula. This is true in 
the present introductory problem, but in a more general problem setting, we 
have some differential equation where SymPy perhaps can help with finding 
an exact solution only in a special case. We can use this special-case solution 
to control our reasoning about scales in the more general setting. 


Theory. The challenges in our scaling is to find the right ue and te scales. 
From (2.8) we see that if u’ — 0 as t + co, u approaches the constant value 
b/a. It can be convenient to let the scaled u > 1 as we approach the du/dt =0 
state. This idea points to choosing 


b o ) p 
TER 1 2.10 
a 9 (2 V (PN 


On the sign of the scaled velocity 


A little note on the sign of ue is necessary here. With op < o, the buoy- 
ancy force upwards wins over the gravity force downwards, and the body 


will move upwards. In this case, the terminal velocity ue > 0. When 
0» > 0, we get a motion downwards, and ue < 0. The corresponding u 
is then also negative, but the scaled velocity u/uc, becomes positive. 


Inserting u = uct = bu/a and t = tet in (2.8) leads to 


du te i a 
= = —t,au4 ; =I-. 
Ji tau Fie u(0) 5 
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We want the scales such that du/dt and ŭ are about unity. To balance the 
size of ù and du/dt we must therefore choose te = 1/a, resulting in the scaled 
ODE problem 


dia 
= ath, u(0) = B, (2.11) 
where 8 is a dimensionless number, 
I 
paiar (2.12) 
Ue b 


reflecting the ratio of the initial velocity and the terminal (t > o0) velocity 

b/a. Scaled equations normally end up with one or more dimensionless param- 

eters, such as 8 here, containing ratios of physical effects in the model. Many 

more examples on dimensionless parameters will appear in later sections. 
The analytical solution of the scaled model (2.11) reads 


elt) =e™ (è —14+ 8) =14+(8-le*. (2.13) 


The result (2.11) with the solution (2.13) is actually astonishing if a and 
b are as in (2.9): the six parameters d, 11, 0b, 0, V, and g are conjured to one: 


dul = 
enV g \ eb 


which is an enormous simplification of the problem if our aim is to investigate 
how u varies with the physical input parameters in the model. In particular, 
if the motion starts from rest, 6 = 0, and there are no physical parameters 
in the scaled model! We can then perform a single simulation and recover all 
physical cases by the unscaling procedure. More precisely, having computed 
u(t) from (2.11), we can use 


u(t) = lafat), (2.14) 


to scale back to the original problem again. We observe that (2.11) can utilize 
a solver for (2.8) by setting a = 1, b = 1, and I = 8. Given some implementa- 
tion of a solver for (2.8), say solver(I, a, b, T, dt, theta), the scaled 
model is run by solver (beta, 1, 1, T, dt, theta). 


Software. We may develop a solver for the scaled problem that uses joblib 
to cache solutions with the same (6, At, and T. For now we fix 6 = 0.5. 
The module decay_vc.py (see Section 3.1.3 in [3] for details) has a function 
solver(I, a, b, T, dt, theta) for solving u’(t) = —a(t)u(t)+0(t) for t € 
(0,T], u(0) = I, with time step dt. We reuse this function and call it with 
a=b=1and I = to solve the scaled problem: 
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from decay_vc import solver as solver_unscaled 


def solver_scaled(beta, T, dt, theta=0.5): 
Solve u’=-u+1, u(0)=beta for (0,T] 
with step dt and theta method. 
print ’Computing the numerical solution’ 
return solver_unscaled( 
I=beta, a=lambda t: 1, b=lambda t: 1, 
T=T, dt=dt, theta=theta) 


import joblib 
disk_memory = joblib.Memory(cachedir=’ temp’) 
solver_scaled = disk_memory.cache(solver_scaled) 


If we want to plot the physical solution, we need an unscale function, 


def unscale(u_scaled, t_scaled, d, mu, rho, rho_b, V): 
a, b = ab(d, mu, rho, rho_b, V) 
return (b/a)*u_scaled, a*t_scaled 


def ab(d, mu, rho, rho_b, V): 
g = 9.81 
a = 3*pi*d+*mu/(rho_b*V) 
b = g*(rho/rho_b - 1) 
return a, b 


Looking at droplets of water in air, we can fix some of the parameters 
and let the size parameter d be the one for experimentation. The following 
function sets physical parameters, computes 3, runs the solver for the scaled 
problem (joblib detects if it is necessary), and finally plots the scaled curve 


u(t) and the unscaled curve u(t). 


def main(dt=0.075, # Time step, scaled problem 


T=7.5, # Final time, scaled problem 
d=0.001, # Diameter (unscaled problem) 

I=0, # Initial velocity (unscaled problem) 
DE: 


# Set parameters, solve and plot 

rho = 0.00129E+3 # air 

rho_b = 1E+3 # density of water 

mu = 0.001 # viscosity of water 

# Asumme we have list or similar for d 

if not isinstance(d, (list,tuple,np.ndarray)): 


aita] 
legends1 = [] 
legends2 = [] 


plt.figure(1) 
plt.figure(2) 
betas = [] # beta values already computed (for plot) 


30 2 Ordinary differential equation models 


for d_ in d: 
V = 4*pi/3*(d_/2.)**3 # volume 
Ely le = a 5 fail, g aeleVe)_loy, VO 
beta = I*a/b 
# Restrict to 3 digits in beta 
beta = abs(round(beta, 3)) 


print ’beta=%.3f’ % beta 
u_scaled, t_scaled = solver_scaled(beta, T, dt) 


# Avoid plotting curves with the same beta value 
if not beta in betas: 

plt.figure(1) 

plt.plot(t_scaled, u_scaled) 

plt.hold(’on’) 

legends1.append(’beta=/g’ % beta) 
betas. append (beta) 


plt.figure(2) 
u, t = unscale@iscailled, t scaled da mu, rho, cho b,, V) 
plt.plot(t, u) 
plt.hold(’on’) 
legends2.append(’d=%g¢ [mm]’ % (d_*1000)) 
plt.figure(1) 
plt.xlabel(’scaled time’); plt.ylabel(’scaled velocity’) 
plt.legend(legends1, loc=’lower right’) 


can be skipped when trying to understand how we work with a scaled model 
to perform the computations. The complete program is found in the file 
falling body.py. 

Since J = 0 implies 6 = 0, we can run different d values without any need 
to recompute u(t) as long as we assume the particle starts from rest. 

From the scaling, we see that ue = b/a ~ d7? and also that te = 1/a ~ d7?, 
so plotting of u(t) with dimensions for various d values will involve significant 
variations in the time and velocity scales. Figure 2.2 has an example with 
d=1,2,3 mm, where we clearly see the different time and velocity scales in 
the figure with unscaled variables. Note that the scaled velocity is positive 
because of the sign of ue (see the box above). 


1.0, 0 


scaled velocity 
u [m/s] 


oe =a — d=1 [mm] 
— d=2 [mm] 
— beta=0 — d=3 [mm] 


1 2 3 4 5 6 7 8 0 20 40 60 80 100 120 140 
scaled time t[s] 


Fig. 2.2 Velocity of falling body: scaled (left) and with dimensions (right). 
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2.1.6 Variable coefficients 


When a prescribed coefficient like a(t) in u’(t) = —a(t)u(t) varies with time 
one usually also performs a scaling of this a, 


ag = US. 
Qe 
where the goal is to have the scaled @ of size unity: |a| < 1. This property 
is obtained by choosing ae as the maximum value of |a(t) — ao] for t € [0,7], 
which is usually a quantity that can be estimated since a(t) is known as a 
function of t. The ag parameter can be chosen as 0 here. (It could be tempting 
to choose ag = min; a(t) so that 0 <a < 1, but then there is at least one point 
where a = 0 and the differential equation collapses to u’ = 0.) 

As an example, imagine a decaying cell culture where we at time tı change 
the environment (typically the nutrition) such that the death rate increases 
by a factor 5. Mathematically, a(t) = d for t < tı and a(t) = 5d for t > tı. The 
model reads u’ = —a(t)u, u(0) = I. 

The a(t) function is scaled by letting the characteristic size be ac = d and 
ag = 0: 


Ome Say 


5, t> tı/te 
The scaled equation becomes 


uc du ee = 
a =acA(t)uct, ucu(0) =I. 


The natural choice of ue is J. The characteristic time, previously taken as 
te = 1/a, can now be chosen as te = t1 or te = 1/d. With te = 1/d we get 

u(t)=—au, u(0)=1, a= { = (2.15) 
where 


y=tid 


is a dimensionless number in the problem. With te = t1, we get 


The dimensionless parameter y is now in the equation rather than in the 
definition of a. Both problems involve y, which is the ratio between the time 
when the environmental change happens and the typical time for the decay 


(1/d). 
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A computation with the scaled model (2.15) and the original model with 
dimensions appears in Figure 2.3. 


— gamma=0.833 — d=0.01 [1/s], t_1=100.00s 


08 0.8 


scaled velocity 
u 


0.0 0.0 
0.0 0.5 1.0 15 2.0 2.5 3.0 3.5 4.0 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0,035 


scaled time t[s] 


Fig. 2.3 Exponential decay with jump: scaled model (left) and unscaled model (right). 


2.1.7 Scaling a cooling problem with constant 
temperature in the surroundings 


The heat exchange between a body at temperature T(t) and the surroundings 
at constant temperature Ts can be modeled by Newton’s law of cooling: 


T'(t)=—k(T—-T;), T(0)= To, (2.16) 
where k is a prescribed heat transfer coefficient. 


Exact solution. An analytical solution is always handy to have as a control 
of the choice of scales. The solution of (2.16) is by standard methods for 
ODEs found to be T(t) = Ts + (To —Ts)e7™. 


Scaling. Physically, we expect the temperature to start at To and then to 
move toward the temperature of the surroundings (Ts). We therefore expect 
that T lies between To and T,. This is mathematically demonstrated by the 
analytical solution as well. A proper scaling is therefore to scale and translate 
T according to 


(2.17) 


Now, 0<T <1. g 
Scaling time by t = t/te and inserting T = To + (Ts —To)T and t = tet in 
the problem (2.16) gives 
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dT = = 
— = —tek(T —1 T(0)=0. 
= = tek(—1), TO) 
A natural choice, as argued in other exponential decay problems, is to choose 
tek = 1, which leaves us with the scaled problem 
dT = = 
—=-(T-1), T(0)=0. (2.18) 
dt 
No physical parameter enters this problem! Our scaling implies that T starts 
at 0 and approaches 1 as t — 00, also in the case Ts < To. The physical 
temperature is always recovered as 


T(t) = To + (Ts —T)T (kt). (2.19) 


Software. An implementation for (2.16) works for (2.18) by setting k = 1, 
T;=1, and Ty =0. 


Alternative scaling. An alternative temperature scaling is to choose 


(2.20) 


Now T = 1 initially and approaches zero as t > 00. The resulting scaled ODE 
problem then becomes 


T(0) =1,. (2.21) 


with solution T = e7*. 


2.1.8 Scaling a cooling problem with time-dependent 
surroundings 


Let us apply the model (2.16) to the case when the surrounding temperature 
varies in time. Say we have an oscillating temperature environment according 
to 


T(t) = Tm + asin(wt), (2.22) 


where Tm is the mean temperature in the surroundings, a is the amplitude 
of the variations around Tm, and 27/w is the period of the temperature 
oscillations. 


Exact solution. Also in this relatively simple problem it is possible to solve 
the differential equation problem analytically. Such a solution may be a good 
help to see what the scales are, and especially to control other forms for 
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reasoning about the scales. Using the method of integrating factors for the 
original differential equation, we have 


t 
T(t) = Tye "t +e | eT, (r)dr. 
0 
With T;,(t) = Tm + asin(wt) we can use SymPy to help us with integrations 
(note that we use w for w in the computer code): 
>>> from sympy import * 


>>> t, k, T_m, a, w = symbols(’t k T_m a w’, real=True, positive=True) 
>>> T_s = Tm + a*sin(w*t) 


>>> I = exp(k*t)*T_s 

>>> I = integrate(I, (t, 0, t)) 
>>> Q = k*exp(-k+t) *I 

>>> Q = simplify (expand (Q)) 


>>> print Q 

(-T_m*k**2 - T_m*w**2 + a*k*w + 

(T_m*k**2 + T_m*w**2 + axk**2*sin(t*w) - 

axk*w*cos (t*w) ) *exp(k*t) )*exp(-k+*t)/((k**2 + w**2)) 


Reordering the result, we get 


T(t) = Toe"! FTm(l e**) + (k?4 w?) t (akwe™ "t + aksin(wt) —akwcos(wt)). 


Scaling. The scaling (2.17) brings in a time-dependent characteristic tem- 
perature scale T; — To. Let us start with a fixed scale, where we take the 
characteristic temperature variation to be Tm — To: 
Fa T -To l 
Tm = To 

We realize by physical reasoning that T sets out at Tọ, but with time, it 
will oscillate around Tm. (This reasoning can be controlled by looking at the 
exact solution we produced above.) The typical average temperature span is 
therefore |Tm — Tol, unless a is much larger than |Tm —To| or To is very close 
to Tim. 

We get from the differential equation, with te = 1/k as in the former case, 


dT = 
k(Tm — To) P k((Tm — To)T + To — Tm — asin(wt), 
resulting in 
ar — l 7 
o —T+1+asin(6t), T(0)=0, (2.23) 


where we have two dimensionless numbers: 
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The a quantity measures the ratio of temperatures: amplitude of oscillations 
versus distance from initial temperature to the average temperature for large 
times. The 8 number is the ratio of the two time scales: the frequency of the 
oscillations in T, and the inverse e-folding time of the heat transfer. For clear 
interpretation of 8 we may introduce the period P = 27/w of the oscillations 
in T, and the e-folding time e = 1/k. Then 8 = 27e/P and measures the 
e-folding time versus the period. 


Remark 


The original problem features five physical parameters: k, To, Tm, a, 
and w, but only two dimensionless numbers appear in the scaled model 
(2.23). In fact, this is an example where application of the Pi theorem 


(see Section 1.1.3) falls short. Since, only time and temperature are 
involved as unit types, the theorem predicts that the five parameters 
yields three dimensionless numbers, not two. Scaling of the differential 
equations, on the other hand, shows us that the two parameters Tm and 
To affect the nature of the problem only through their difference. 


Software. Implementations of the unscaled problem (2.16) can be reused for 
the scaled model by setting k = 1, To =0, and T,(t) = 1 +asin(bt) (Tn, =1, 
a=a,w =). The file osc_cooling.py contains solvers for the problem with 
dimensions and for the scaled problem. The figure below shows three cases 
of @ values: small, medium, and large. 


1.2 


For the small 8 value, the oscillations in the surrounding temperature are 
slow enough compared to k for the heating and cooling process to follow the 
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surrounding temperature, with a small time lag. For larger 8, the heating 
and cooling require more time, and the oscillations get smaller. 


Discussion of the time scale. There are two time variations of importance 
in the present problem: heat is transferred to the surroundings at a rate k, 
and the surroundings have a temperature variation with a period that goes 
like 1/w. (We can, when we are so lucky that we have an analytical solution 
at hand, inspect this solution to see that k impacts the problem through a 
decay factor e~**, and w impacts the problem through oscillations sin(wt).) 
The k parameter related to temperature decay points to a time scale te =1/k, 
while the temperature oscillations of the surroundings point to a time scale 
te =1/w. Which one should be chosen? 

Bringing the temperature from Tp to the level of the surroundings, Tm, 
goes like e~**, so in this process te = 1/k is the characteristic time. Thereafter, 
the body’s temperature just responds to the oscillations and the sin(wt) (and 
cos(wt)) term dominates. For these large times, te = 1/w is the appropriate 
time scale. Choosing te = 1/w results in 

dT 
dt 

Let us illustrate another, less effective, scaling. The temperature scale 
in (2.17) looks natural, so we apply this choice of scale. The characteristic 


temperature To — T, now involves a time-dependent term T;(t). The mathe- 
matical steps become a bit more technically involved: 


—p-\(T—(1+asin(#))), T(0)=0. (2.24) 


aT _ dTs 5) dT dt 
dt a dt T T (Ts To) 


With t=t/t. = kt we get from the differential equation 


dT, = dT p 
=P A(T hy) k= -kT -IT-T 
dt + (Ts — To) a ( )(Ts — To), 


which after dividing by k(Ts — To) results in 


dT - dTs T 
dt (PE dt k(T,—Tp’ 
or 
T = t/k = 
dT _ (T-1) a ) 
dt k(Tm + asin(wt/k) — To) 


The last term is complicated and becomes more tractable if we factor out 
dimensionless numbers. To this end, we scale T; by (e.g.) Tm, which means 
to factor out Tm in the denominator. We are then left with 
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dT = cos(bt) > 

a(i = T 2.25 
Oe TS eee saa (225) 
where a, 3, and y are dimensionless numbers characterizing the relative im- 
portance of parameters in the problem: 


a=a/Tm, B=w/k, y=To/Tm. (2.26) 


We notice that (2.25) is not a special case of the original problem (2.16). 
Furthermore, the original five parameters k, Tm, a, w, and To are reduced 
to three dimensionless parameters. We conclude that this scaling is inferior, 
because using the temperature scale To — Tm enables reuse of the software 
for the unscaled problem and only two dimensionless parameters appear in 
the scaled model. 

Let us briefly mention another possible temperature scaling: T = T/Tm, 
motivated by the fact that as t > 00, T will oscillate around Tm, so this T 
will oscillate around unity. We get the dimensionless ODE 

GT ek E ESE 

dt 
with a new dimensionless parameter 6 = a/Tm. However, the initial condi- 
tion becomes T(0) = To/Tm, and the ratio To/Tm is a third dimensionless 
parameter, so this scaling is also inferior to the one above with only two 
parameters. 


2.1.9 Scaling a nonlinear ODE 


Exponential growth models, u’ = au, are not realistic in environments with 
limited resources. However, by letting a depend on u, the effect of limited 
resources can well be captured by such a simple differential equation model: 


uw =a(u)ju, u(0)=T. (2.27) 


If the maximum value of u is denoted by M, we have that a(M)=0. A 
simple choice fulfilling this requirement is a(w) = e(1—u/M). The parameter 
o can be interpreted as the initial exponential growth rate if we assume that 
I/M «1, since at t= 0 the model then approximates u’ = ou. 

The choice a(u) = e(1—u/M) is known as the logistic model for population 
growth: 


u' = ou(1—u/M), u(0)=T. (2.28) 


A more complicated choice of a may be a(u) = e(1—u/M)? for some exponent 
p (this function also fulfills a( M) = 0 and a % o for t=0). 
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Scaling. Let us scale (2.27) with a(u) = o(1 — u/M)P. The natural scale 
for u is M (ue = M), since we know that 0 < u < M, and this makes the 
dimensionless i = u/M € (0,1]. The function a(u) is typically varying between 
0 and ọ, so it can be scaled as 


aa) =“) a- Bp =a. 


Time is scaled as t= t /te for some suitable characteristic time te. Inserted in 
(2.27), we get 


uc du Si oe = 
d = oūucŭ, ucU0)=I, 
resulting in 
du Cip 4 I 
T =to(l1— u) u, u(0)= i 


A natural choice is te = 1/0 as in other exponential growth models since 
it leads to the term on the right-hand side to be about unity, just as the 
left-hand side. (If the scaling is correct, u and its derivatives are of order 
unity, so the coefficients must also be of order unity.) Introducing also the 
dimensionless parameter 


a= M’ 
measuring the fraction of the initial population compared to the maximum 
one, we get the dimensionless model 


du 

dt 
Here, we have two dimensionless parameters: a and p. A classical logistic 
model with p = 1 has only one dimensionless variable. 


=(1-a)?a, u(0)=a. (2.29) 


Alternative scaling. We could try another scaling of u where we also trans- 
late u: 


oes u—-L 
u= M 
This choice of u& results in 
du om) 4. 
ee (l-a—u)Pu, u(0)=0. (2.30) 


The essential difference between (2.29) and (2.30) is that wu € [a,1] in the 
former and ù € [0,1 — a] in the latter. Both models involve the dimensionless 
numbers a and p. An advantage of (2.29) is that software for the unscaled 
model can easily be used for the scaled model by choosing J =a, M = 1, and 
o=1. 
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2.1.10 SIR ODE system for spreading of diseases 


The field of epidemiology frequently applies ODE systems to describe the 
spreading of diseases, such as smallpox, measles, plague, ordinary flu, swine 
flu, and HIV. Different models include different effects, which are reflected in 
dimensionless numbers. Most of the effects are modeled as exponential decay 
or growth of the dependent variables. 

The simplest model has three categories of people: susceptibles (S) who can 
get the disease, infectious (I) who are infected and may infect susceptibles, 
and recovered (R) who have recovered from the disease and gained immunity. 
We introduce S(t), I(t), and R(t) as the number of people in the categories 
S, I, and R, respectively. The model, naturally known as the SIR model!, can 
be expressed as a system of three ODEs: 


> = —$SI, (2.31) 
dI 

= =pSI-vI, (2.32) 
=v], (2.33) 


where 6 and v are empirical constants. The average time for recovering from 
the disease can be shown to be v~!, but 8 is much harder to estimate, so 
working with a scaled model where ĝ is “scaled away” is advantageous. 


Scaling. It is natural to scale S, I, and R by, e.g., S(0): 
S - I = R 
I= R= —. 
S(0)’ S(0)’ S(0) 


Introducing t = t/t, we arrive at the equations 


S= 


dS A 
se ip 

FE = tes O5 

di Sz - 
nae Pat, 
zz = teS(0)88 7 
dR = 

a = tl, 

am as 


with initial conditions $(0) = 1, I(0) = Io/S(0) =a, and R(0) = R(0)/S(0). 
Normally, R(0) = 0. 

Taking te = 1/v, corresponding to a time unit equal to the time it takes 
to recover from the disease, we end up with the scaled model 


lhttps://en.wikipedia.org/wiki/Epidemic_model 
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dS z. 

— = —RoSI 2.34 
dt 0 , ( ) 
dI z= = 

— = I-I 2: 

== Rob] -], (2.35) 
dR - 

~= 2. 
RL, (2.36) 


with $(0) = 1, 7(0) =a, R(0) = 0, and Ro as the dimensionless number 
S(0)6 


v 
We see from (2.35) that to make the disease spreading, dI /dt > 0, and there- 
fore RoS(0)—1>0 or Ro > 1 since S(0) = 1. Therefore, Ro reflects the 
disease’s ability to spread and is consequently an important dimensionless 
quantity, known as the basic reproduction number?. This number reflects the 
number of infected people caused by one infectious individual during the time 
period of the disease. 

Looking at (2.32), we see that to increase J initially, we must have dI /dt > 0 
at t= 0, which implies 8I (0)S(0)—vI(0) > 0, i.e., Ro > 1. 


Ro= (2.37) 


Software. Any implementation of the SIR model with dimensions can be 
reused for the scaled model by setting 8 = Ro, v = 1, $(0) =1—a, and 
I(0) = a. Below is a plot with two cases: Ro = 2 and Rọ = 5, both with 
a= 0.02. 


RO=2, alpha=0.02 


RO=5, alpha=0.02 


-uv 


Alternative scaling. Adding (2.31)-(2.33) shows that 


dS dI dR 
nru a | => S+I+R=const=N, 


*https://en.wikipedia. org/wiki/Basic_reproduction_number 
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where N is the size of the population. We can therefore scale S, J, and R by 
the total population N = $(0)+J(0)+ R(0): 
= S&S >- I - R 
S= N’ p= N’ R= N 
With the same time scale, one gets the system (2.34)-(2.36), but with Ro 
replaced by the dimensionless number: 
~ N 
ee (2.38) 
v 
The initial conditions become $(0) = 1 — a, I(0) = a, and R(0) = Oy. > 
For the disease to spread at t= 0, we must have RpS(0) > 1, but RoS(0) = 
N£6/v-S(0)/N = Ro, so the criterion is still Ro > 1. Since Ro is a more famous 
number than Ro, we can write the ODEs with Ro/S(0) = Ro/(1—a) instead 
of Ro. 
Choosing te to make the SJ terms balance the time derivatives, te = 
(NB)~1, moves Ro (or Ro if we scale S, I, and R by S(0)) to the J terms: 


2.1.11 SIRV model with finite immunity 


A common extension of the SIR model involves finite immunity: after some 
period of time, recovered individuals lose their immunity and become suscep- 
tibles again. This is modeled as a leakage — uR from the R to the S category, 
where u`! is the average time it takes to lose immunity. Vaccination is an- 
other extension: a fraction pS' is removed from the S category by successful 
vaccination and brought to a new category V (the vaccinated). The ODE 
model reads 
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dS 


Fp = ASE pS + wR, (2.39) 
H = 6SI—v1, (2.40) 
a ae (2.41 
BV as (2.42) 


Using te = 1/v and scaling the unknowns by $(0), we arrive at the dimen- 
sionless model 


dS =u, z 

> = —RoSI — ôS +4R, (2.43) 
dI fee e 

— = I—II 2.44 
dt RoS oy ( ) 
dR - >- 

— =] -yR 2.45 
d yR, (2.45) 
dV = 

Lan 2.4 
= = 85, (2.46) 


with two new dimensionless parameters: 


y=", =Ë. 
V V 


The quantity p™ t can be interpreted as the average time it takes to vaccinate 


a susceptible successfully. Writing y = v~!/u-! and 6=v~!/p~! gives the 
interpretation that y is the ratio of the average time to recover and the 
average time to lose immunity, while 6 is the ratio of the average time to 
recover and the average time to successfully vaccinate a susceptible. 

The plot in Figure 2.4 has y = 0.05, i.e., loss of immunity takes 20 weeks 
if it takes one week to recover from the disease. The left plot corresponds to 
no vaccination, while the right has 6 = 0.5 for a vaccination campaign that 
lasts from day 7 to day 15. The value 6 = 0.5 reflects that it takes two weeks 
to successfully vaccinate a susceptible, but the effect of vaccination is still 
dramatic. 


2.1.12 Michaelis-Menten kinetics for biochemical 
reactions 


A classical reaction model in biochemistry describes how a substrate S is 
turned into a product P with aid of an enzyme E. S and E react to form 
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Ry =5, a=0.02, y=0.05, 6=0 R, =5, a=0.02, y=0.05, 6=0.5 


<27 u 


“Oo 10 20 30 40 50 60 


Fig. 2.4 Spreading of a disease with loss of immunity (left) and added vaccination 
(right). 


a complex ES in the first stage of the reaction. In the second stage, ES is 
turned into E and P. Introducing the amount of S, E, ES, and P by [S], [E], 
[ES], and [P], the mathematical model can be written as 


aes) = k4[E][S] —k,[ES]—k_[BS], (2.47) 
d[P 

as) = —k,[E][5] + k_[BSI, (2.49) 
a) = —ky[E][S]+k_[BS] + ky [BS]. (2.50) 


The initial conditions are [E£'S](0) = [P](0) = 0, and [S] = So, [E] = Eo. Three 
rate constants are involved: k4, k_, and ky. The above mathematical model 
is known as Michaelis-Menten kinetics®. 

The amount of substance is measured in the unit mole* (mol). From the 
equations we can see that k is measured in s~!mol~!, while k— and ky are 
measured in s71. It is convenient to get rid of the mole unit for the amount 
of a substance. When working with dimensionless quantities, only ratios of 
the rate constants and not their specific values are needed. 


Classical analysis. A common assumption is that the formation of [E'S] is 
very fast and that it quickly reaches an equilibrium state, [ES]! = 0. Equation 
(2.47) then reduces to the algebraic equation 

k4[E][$]— kv [ES]—k_[ES] =0, 
which leads to 


3https://en.wikipedia.org/wiki/Michaelis-Menten_kinetics 
fhttps ://en.wikipedia.org/wiki/Mole_(unit) 
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[E][S] k-t+ky 
= =K, 2.51 
EJ k eni 
where K is the famous Michaelis constant - the equilibrium constant between 
[E][S] and [ES]. 
Another important observation is that the ODE system implies two con- 
servation equations, arising from simply adding the ODEs: 


+4 =0, (2.52) 


dt ae 239) 

from which it follows that 
[ES] + [E] = Eo, (2.54) 
[ES|+[S]+[P] = So. (2.55) 


We can use (2.54) and (2.51) to express [E] by [S]: 


KE 
= Tg 


Now (2.49) can be developed to an equation involving [S] only: 


[E] = Eo- [ES] = Eo - E2 


ASI k [EJS] +k- [ES] 
k_ 
= (-k; + =) BI5) 
k_ KE 
=(-k + SET Tg 
=- y (2.56) 


We see that the parameter K is central. 
From above expression for |E] and (2.54) it now follows 


K Eo Eo[S] 
E\=——,,,_ [ES] = =. 
[E] Kris) | | K+[S] 
If K is comparable to So these indicate 
Eo So 


(E]~ Eo [ES]~ 2, 
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as is used for scaling |E] and Qc, subsequently. Provided we exclude the case 
[S] > K, we may infer that [E] will be of magnitude Eo, while [ES] will be 
of magnitude EgSo/K. 


Dimensionless ODE system. Let us reason how to make the original ODE 
system (2.47)-(2.50) dimensionless. Aiming at [S] and [E] of unit size, two 
obvious dimensionless unknowns are 


go) ge), 
So Eo 
For the other two unknowns we just introduce scales to be determined later: 
pakl oe 
Pe Qe 


With t= t/t. the equations become 


= aon = BS —te(ky +h_)Q, 
af = teks SEQ, 

=a = -toky Eo ES + tek- Z 3, 

= = -tek SoES + te(k— +h) 22. 


Determining scales. Choosing the scales is actually a quite complicated 
matter that requires extensive analysis of the equations to determine the 
characteristics of the solutions. Much literature is written about this, but 
here we shall take a simplistic and pragmatic approach. Besides the Michaelis 
constant K, there is another important parameter, 


because most applications will involve a small e. We shall have K and e in 
mind while choosing scales such that these symbols appear naturally in the 
scaled equations. 

Looking at the equations, we see that the K parameter will appear if 
te~ 1/k4.. However, 1/k4 does not have the dimension [T]~' as required, so 
we need to add a factor with dimension mol. A natural choice is t7 t = k4 Sp 
or t71 = k} Eo. Since often So >> Eo, the former te is a short time scale and 
the latter is a long time scale. If the interest is in the long time scale, we set 


a2 a 
~ ky Eo” 


te 
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The equations then take the form 


dQ So -1 

— = — ES- KE 

dt Qe 3 0 Q 

dP ky Qe > 

dt kyo Pe 

dS an en a 
—-=—ES 4 ; 

dt k4 Eo So Q 

dE -lpo Qe A 

Seo ERKE 
dt : Be” 


The [ES] variable starts and ends at zero, and its maximum value can be 
roughly estimated from the equation for [ES] by setting [ES] = 0, which 
gives 
[E][S] _ EoSo 

ES] = — => ~ —— 

(Bs) = a 00, 
where we have replaced [E][S] by EoSo as an identification of magnitude. 
This magnitude of [ES] at its maximum can be used as the characteristic 
size Qe: 


Eo So 
K 


The equation for P simplifies if we choose P, = Qe. With these assumptions 
one gets 


Qc= 


dQ _ K-(BS 6), 
dt 

dP ke = 

dt kE | 

dS. ge k- Eos 
di ERK” 
“ =-e 1£S+e'9 


where we see that a= 8 +y, so y can be eliminated. Moreover, 
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k— 


= E Eo 


+8, 


implying that a> p. 
We arrive at the final set of scaled differential equations: 


dQ =a = 
AREN (2.57) 
dP = 
dE = BQ, (2.58) 
i = -E8 +(1-—ßa7t)Q, (2.59) 
dE Sa 

ee = —ES+Q. (2.60) 


The initial conditions are S = E = 1 and Q = P=0. 
The five initial parameters (So, Eo, ki, k—, and ky) are reduced to three 
dimensionless constants: 


e ais the dimensionless Michaelis constant, reflecting the ratio of the pro- 
duction of P and E (k, +k_) versus the production of the complex (k+), 
made dimensionless by Eo, 

e ¢ is the initial fraction of enzyme relative to the substrate, 

e 8 measures the relative importance of production of P (ky) versus produc- 
tion of the complex (k), made dimensionless by Eo. 


Observe that software developed for solving (2.47)-(2.50) cannot be reused for 
solving (2.57)-(2.60) since the latter system has a slightly different structure. 


Conservation equations. The counterpart to the conservation equations 
(2.54)-(2.55) is obtained by adding (2.57) and a times (2.60), and adding 
(2.57), (2.58), and a times (2.59): 


«ta 1Q+E=1, (2.61) 
aS+Q+P=a. (2.62) 


The scaled quantities, as well as the original concentrations, must be pos- 
itive variables, and E € [0,1], S € [0,1]. Such checks along with the conserved 
quantities above should be performed at every time step in a simulation. 


Analysis of the scaled system. In the scaled system, we may assume e€ 
small, which from (2.60) gives rise to the simplification «E’ = 0, and thereby 
the relation Q = ES. The conservation equation [ES] + [E] = Ep reads Q-Q + 
EoE = Ep such that E = 1—Q,Q/Ep = 1-—QSp/K = 1 — eta tQ. The 
relation Q = ES then becomes 
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which can be solved for Q: 


The equation (2.59) for S becomes 
dS = 5 
N | A ee 
dt ate 1$ 

This is a more precise analysis than the one leading to (2.56) since we now 

realize that the mathematical assumption for the simplification is € > 0. 

Is (2.63) consistent with (2.56)? It is easy to make algebraic mistakes when 
deriving scaled equations, so it is always wise to carry out consistency checks. 
Introducing dimensions in (2.63) leads to 


(2.63) 


tedS dS BS ky S ky S 


Sodt dt  a+e!5  k}Eo KE +E o5 ky K+S 


and hence with t7! = ki Eo, 


dS kyEyS 


dt K+S’ 


which is (2.56). 

Figure 2.5 shows the impact of e: with a moderately small value (0.1) we 
see that Q ~ 0, which justifies the simplifications performed above. We also 
observe that all the unknowns vary between 0 and about 1, indicating that the 
scaling is successful for the chosen dimensionless numbers. The simulations 
made use of a time step At = 0.01 with a 4th-order Runge-Kutta method, 
using a = 1.5, 6 = 1 (relevant code is in the simulate_biochemical_process 
function in session. py). 


alpha=1.5, beta=1, epsilon=0.9 


alpha=1.5, beta=1, epsilon=0.1 


0.8 


complex || 0.6 
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Fig. 2.5 Simulation of a biochemical process. 
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However, it is of interest to investigate the limit e — 0. Initially, the equa- 
tion for dE/dt reads dE/dt = —e~, which implies a very fast reduction of 
E. Using e = 0.005 and At = 1078, simulation results show that E decays to 
approximately zero at t = 0.03 while S ~ 1 and Q ~ P ~0. This is reasonable 
since with very little enzyme in comparison with the substrate (e — 0) very 
little will happen. 


2.2 Vibration problems 


We shall in this section address a range of different second-order ODEs for 
mechanical vibrations and demonstrate how to reason about the scaling in 
different physical scenarios. 


2.2.1 Undamped vibrations without forcing 


The simplest differential equation model for mechanical vibrations reads 


mu’+ku=0, wu(0)=T, u/(0)=V, (2.64) 


where unknown u(t) measures the displacement of the body, This is a common 
model for a vibrating body with mass m attached to a linear spring with 
spring constant k (and force —ku). Figure 2.6 shows a typical mechanical 
sketch of such a system: some mass can move horizontally without friction 
and is connected to a spring that exerts a force —ku on the body. 


b>) 


Fig. 2.6 Oscillating body attached to a spring. 


The first technical steps of scaling. The problem (2.64) has one inde- 
pendent variable t and one dependent variable u. We introduce dimensionless 
versions of these variables: 
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_ u - t 
u=—, t=, 
Uc te 
where uc and te are characteristic values of u and t. Inserted in (2.64), we 
get 


Ue au = = Ue du 
resulting in 
da tk Fang) = Me 
á 2 = — T = . 2 
e naa a0)=—, w= — (2.65) 


What is an appropriate displacement scale uc? The initial condition u(0) = 
I is a candidate, i.e., ue = I. But how to choose the time scale? Making the 
coefficient in front of the ù unity, such that both terms balance and are of 
size unity, is a candidate. 


The exact solution. To better see what the proper scales of u and t are, 
we can look into the analytical solution of this problem. Although the exact 
solution of (2.64) is quite straightforward to calculate by hand, we take the 
opportunity to make use of SymPy to find u(t). The use of SymPy can later 
be generalized to vibration ODEs that are harder to solve by hand. 

SymPy requires all mathematical symbols to be explicitly created: 


from sympy import * 

u = symbols(’u’, cls=Function) 

w = symbols(’w’, real=True, positive=True) 

I, V, C1, C2 = symbols(’I V C1 C2’, real=True) 


To specify the ODE to be solved, we can make a Python function returning 
all the terms in the ODE: 


# Define differential equation: u’? + w**2*u = 0 
def ode(u): 
return diff(u, t, t) + w**2*u 


diffeq = ode(u(t)) 


The diffeq variable, defining the ODE, can be passed to the SymPy function 
dsolve to find the symbolic solution of the ODE: 


s = dsolve(diffeq, u(t)) 

# s is an u(t) == expression (Eq obj.), s.rhs grabs the expression 
u_sol = s.rhs 

print u_sol 


The solution that gets printed is Cl*sin(t*w) + C2*cos(t*w), indicating 
that there are two integration constants C1 and C2 to be determined by the 
initial conditions. The result of applying these conditions is a 2 x 2 linear 
system of algebraic equations that SymPy can solve by the solve function. 
The code goes as follows: 
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# The solution u_sol contains integration constants C1 and C2 
# but these are not symbols, substitute them by symbols 
n oL = m EoLA EC 1) .subs E a E 


# Determine C1 and C2 from the initial conditions 

ic = [u_sol.subs(t, 0) - I, u_sol.diff(t).subs(t, 0) - VJ] 
print ic # 2x2 algebraic system for C1 and C2 

s = solve(ic, [C1, C2]) 

# s is now a dictionary: {C2: I, C1: V/w} 

# substitute solution back in u_sol 

u_sol = u_sol.subs(Ci, s[C1i]).subs(C2, s[C2]) 

print u_sol 


The u_sol variable is now I*cos(t*w) + V*sin(t*w)/w. Since symbolic 
software is far from bug-free and can give wrong results, we should always 
check the answer. Here, we insert the solution in the ODE to see if the result 
is zero, and we insert the solution in the initial conditions to see that these 
are fulfilled: 


# Check that the solution fulfills the ODE and init.cond. 
print simplify(ode(u_sol)), 
print u_sol.subs(t, 0) - I, diff(u_sol, t).subs(t, 0) - V 


There will be many more examples on using SymPy to find exact solutions 
of differential equation problems. 
The solution of the ODE in mathematical notation is 


V. k 
u(t) = I cos(wt) + F sin(wt), w= a 


More insight arises from rewriting such an expression in the form Acos(wt — 


@): 
u(t) = 4/1? + Y cos(wt 9) o=tan!(V/(w1)). 


Now we see that the u corresponds to cosine oscillations with a frequency 
shift ¢ and amplitude yI? + (V/w)?. 


Discussion of the displacement scale. The amplitude of u is y1? + V? /w?, 
and this expression is obviously a candidate for ue. However, the simpler 
choice ue = max(I,V/w) is also relevant and more attractive than the square 
root expression (but potentially a factor 1.4 wrong compared to the exact 
amplitude). It is not very important to have |u| <1, the point is to avoid |u| 
very small or large. 


Discussion of the time scale. What is an appropriate time scale? Looking 
at (2.65) and arguing that u” and u both should be around unity in size, the 
coefficient t?k/m must equal unity, implying that te = ,/m/k. Also from the 
analytical solution we see that the solution goes like the sine or cosine of wt, 
so 1/w = \/m/k can be a characteristic time scale. Likewise, one period of the 
oscillations, P = 27/w, can be the characteristic time, leading to te = 2r /w. 
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The dimensionless solution. With uc = I and te = ym/k we get the 
scaled model 


—, +u=0, wu(0)=1, w(0)=a, (2.66) 


where a is a dimensionless parameter: 


gavom 
IV 


Note that in case V = 0, we have “scaled away” all physical parameters. The 
universal solution without physical parameters is then u(t) = cost. 
The unscaled solution is recovered as 


u(t) = Iu(./k/mt). (2.67) 
This expressions shows that the scaling is simply a matter of stretching or 
shrinking the axes. 


Alternative displacement scale. Using ue = V/w, the equation is not 
changed, but the initial conditions become 


: I Iw Ifk 3 , 
u(0) ils V V a a ’ u (0) 


With uc = V/w and one period as time scale, te = 27,/m/k, we get the 
alternative model 


Pu 25 = -1 = 
qe i u=0, u(0)=a7, u (0)=2r. (2.68) 


The unscaled solution is in this case recovered by 


u(t) = vy Taer yt. (2.69) 


About frequency and dimensions. The solution goes like coswt, where 
w= /m/k must have dimension 1/s. Actually, w has dimension radians per 
second: rad/s. A radian is dimensionless since it is arc (length) divided by 
radius (length), but still regarded as a unit. The period P of vibrations is a 
more intuitive quantity than the frequency w. The relation between P and 
w is P = 2r/w. The number of oscillation cycles per period, f, is a more 
intuitive measurement of frequency and also known as frequency. Therefore, 
to be precise, w should be named angular frequency. The relation between f 
and T is f =1/T, so f = 2rw and measured in Hz (1/s), which is the unit 
for counts per unit time. 
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2.2.2 Undamped vibrations with constant forcing 


For vertical vibrations in the gravity field, the model (2.64) must also take 
the gravity force —mg into account: 


mu” +ku=-—mg. 


How does the new term —mg influence the scaling? We observe that if there 
is no movement of the body, u” = 0, and the spring elongation matches the 
gravity force: ku = —mg, leading to a steady displacement u = —mg/k. We 
can then have oscillations around this equilibrium point. A natural scaling 
for u is therefore 


u—(—mg/k)  uk+mg 
Ue ~~ kue ` 


u = 


The scaled differential equation with the same time scale as before reads 


du ae t2 t 
= uU = 
dt2 uc g Ue 9» 
leading to 
u 
—+u=0. 
dt? 


The initial conditions u(0) = J and u’(0) = V become, with uc = J, 


T mg du mV 
C gay ks cag) VET 


We see that the oscillations around the equilibrium point in the gravity field 
are identical to the horizontal oscillations without gravity, except for an offset 
mg/(kI) in the displacement. 


2.2.3 Undamped vibrations with time-dependent forcing 


Now we add a transient forcing term F(t) to the model (2.64): 
mu” +ku= F(t), u(0)=I, w(0)=V. (2.70) 
Take the forcing to be oscillating: 
F(t) = Acos(yt). 


The technical steps of the scaling are still the same, with the intermediate 
result 
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u tk t2 I Vt 
“tng e Acos(ytt), a0) =>, w(0)= —*. (2.71) 
dt? m Muc Ue Ue 


What are typical displacement and time scales? This is not so obvious without 
knowing the details of the solution, because there are three parameters (J, V, 
and A) that influence the magnitude of u. Moreover, there are two time scales, 
one for the free vibrations of the systems and one for the forced vibrations 
F(t). 


Investigating scales via analytical solutions. As we have seen already 
several times, having access to an exact solution is very fortunate as it allows 
us to directly examine the scales. Also in the present problem it is possible to 
derive an exact solution. We continue the SymPy session from the previous 
section and perform much of the same steps. Note that we use w for w = y k/m 
in the computer code (to obtain a more direct visual counterpart to w). 
SymPy may get confused when coefficients in differential equations contain 
several symbols. We therefore rewrite the equation with at most one symbol 
in each coefficient (i.e., symbolic software is in general more successful when 
applied to scaled differential equations than the unscaled counterparts, but 
right now our task is to solve the unscaled version). The amplitude A/m in 
the forcing term is of this reason replaced by the symbol A1. 


A, Al, m, psi = symbols (°A A1 m psi’, positive=True, real=True) 
def ode(u): 
return diff (u, t, t) + w**2*u - Al*cos(psi*t) 


diffeq = ode(u(t)) 
u_sol = dsolve(diffeq, u(t)) 
u soll S u soll rhe 


# Determine the constants C1 and C2 in u_sol 

# (first substitute our own declared C1 and C2 symbols, 

# then use the initial conditions) 

Wesolw— ussol.subsi(7Gi 4, (C1) SubS (7622) CE2) 

eqs = [u_sol.subs(t, 0) - I, u_sol.diff(t).subs(t, 0) - V] 
s = solve(eqs, [C1, C2]) 

u sol = u_sol.subs(C1, s[Ci]).subs(C2, s[C2]) 


# Check that the solution fulfills the equation and init.cond. 
print simplify (ode(u_sol)) 

print simplify(u_sol.subs(t, 0) - I) 

print simplify(diff(u_sol, t).subs(t, 0) - V) 


u_sol = simplify(expand(u_sol.subs(A1, A/m))) 
print u_sol 


The output from the last line is 


A/m*cos(psi*t)/(-psi**2 + w**2) + V*sin(t*w)/w + 
(A/m + I*psix*2 - I*w**2)*cos(t*w)/(psix*2 - w**2) 


With a bit of rewrite this expression becomes 
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u(t) = = cos( it) + Y sin(wt) (2 | 1) cos(wt). 
Obviously, this expression is only meaningful for Y 4 w. The case Y = w gives 
an infinite amplitude in this model, a phenomenon known as resonance. The 
amplitude becomes finite when damping is included, see Section 2.2.4. 

When the system starts from rest, I = V = 0, and the forcing is the only 
driving mechanism, we can simplify: 


A A 
u(t) = Py eT ala + UE?) cos(wt) 
= egy Oost) —cos(wt)). 


To gain more insight, cos(¢t) — cos(wt) can be rewritten in terms of the mean 
frequency (Yy +w)/2 and the difference in frequency (Y — w)/2: 


uoi (PHM), em 


showing that there is a signal with frequency (Y +w)/2 whose amplitude has 
a (much) slower frequency (4Y —w)/2. Figure 2.7 shows an example on such 
a signal. 


Fig. 2.7 Signal with frequency 3.1 and envelope frequency 0.2. 
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The displacement and time scales. A characteristic displacement can 
in the latter special case be taken as ue = A/(m(w? — 4?)). This is also a 
relevant choice in the more general case J # 0,V #0, unless J or V is so 
large that it dominates over the amplitude caused by the forcing. With ue = 
A/(m(w? —?)) we also have three special cases: w <1, w > y, and Y ~ w. 
In the latter case we need uc = A/(m(w? — 4?)) if we want |u| < 1. When w 
and w are significantly different, we may choose one of them and neglect the 
smaller. Choosing w means uc = A/k, which is the relevant scale if w >> a. In 
the opposite case, w & Y, Ue = A/(mr?). 

The time scale is dominated by the fastest oscillations, which are of fre- 

quency w or w when these are close and the largest of them when they are 
distant. In any case, we set te = 1/max(@,w). 
Finding the displacement scale from the differential equation. Going 
back to (2.71), we may demand that all the three terms in the differential 
equation are of size unity. This leads to te = \/m/k and ue = At? /m = A/k. 
The formula for ue is a kind of measure of the ratio of the forcing and the 
spring force (the dimensionless number A/(ku¢) would be this ratio). 

Looking at (2.72), we see that if Y <w, the amplitude can be approxi- 

mated by A/(mw?) = A/k, showing that the scale uc = A/k is relevant for an 
excitation frequency w that is small compared to the free vibration frequency 
w. 
Scaling with free vibrations as time scale. The next step is to work out 
the dimensionless ODE for the chosen scales. We first select the time scale 
based on the free oscillations with frequency w, i.e., te = 1/w. Inserting the 
expression in (2.71) results in 


—, +u=ycos(st), u(0)=a, u (0)= £. (2.73) 


I 

a= —, (2.74) 
Vte V 

p e (2.75) 
Uc Wue 
tA A 

y ae ae (2.76) 
te Y% 

pae (2.77) 


We remark that the choice of ue has so far not been made. Several different 
cases will be considered below, and we will see that certain choices reduce 
the number of independent dimensionless variables to three. 

The four dimensionless variables above have interpretations as ratios of 
physical effects: 
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e a: ratio of the initial displacement and the characteristic response uc, 

e £: ratio of the initial velocity and the typical velocity measure u,/te, 

e ~: ratio of the forcing A and the mass times acceleration mu,/t? or the 
ratio of the forcing and the spring force kue 

e 6: ratio of the frequencies or the time scales of the forcing and the free 
vibrations. 


Software. Any solver for (2.71) can be used for (2.73). More details are 
provided at the end of Section 2.2.4. 


Choice of ue close to resonance. Now we shall discuss various choices of 
Uc. Close to resonance, when Y ~w, we may set uc = A/(m(w? — y?)). The 
dimensionless numbers become in this case 


etn i 2 
rae ASE 
V VVkm 2 
B= Wue za A (1 ô J 
A 
=—=1-0° 
Y Bitte , 
poe 
w 
With Y = 0.99w, 6 = 0.99, V = 0, a = y = 1 — ô? = 0.02, we have the problem 
ddu ae 2 oi = 
TE +u = 0.02cos(0.99Ł), wu(0) = 0.02, u’(0) =0 


This is a problem with a very small initial condition and a very small forcing, 
but the state close to resonance brings the amplitude up to about unity, see 
the result of numerical simulations with ô = 0.99 in Figure 2.8. Neglecting a, 
the solution is given by (2.72), which here means A = 1—6?, m=1, w=1, 
p = ô: 
u(t) = 2sin(—0.005t) sin(0.995¢) . 

Note that this is a problem which demands very high accuracy in the numer- 
ical calculations. Using 20 time steps per period gives a significant angular 
frequency error and an amplitude of about 1.4. We used 160 steps per period 
for the results in Figure 2.8. 


Unit size of all terms in the ODE. Using the displacement scale uc = A/k 
leads to (2.73) with 
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dt=0.0392699 


0 200 400 600 800 1000 


Fig. 2.8 Forced undamped vibrations close to resonance. 


1I I 
=u ATR’ 
V Vk 
P= ue Aw 
A 
= = 1 
Y kuc , 
m 
W 


Simulating a case with ô = 0.5, a = 1, and 6 =0 gives the oscillations in 
Figure 2.9, which is a case away from resonance, and the amplitude is about 
unity. However, choosing ô = 0.99 (close to resonance) results in a figure 
similar to Figure 2.8, except that the amplitude is about 10? because of the 
moderate size of ue. The present scaling is therefore most suitable away from 
resonance, and when the terms containing coswt and sinwt are important 
(e.g., w > Ww). 

Choice of ue when w >w. Finally, we may look at the case where ~ >> w 
such that ue = A/ (my?) is a relevant scale (i.e., omitting w? compared to y? 
in the denominator), but in this case we should use te = 1/w since the force 
varies much faster than the free vibrations of the system. This choice of te 
changes the scaled ODE to 
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z0 dt=0.0392699 


50 100 150 


ae +6-7u=ycos(t), u(0)=a, w(0)=8, (2.78) 
where 


SE 
Ue Ajk ’ 
é Vk 
pa V My 
Ue A 
2 
_ tA =1, 
Mue 
te yp 
6= r 
yi w 


In the regime Y >w, ô> 1, thus making a and £8 large. However, if a and/or 
is large, the initial condition dominates over the forcing, and will also 
dominate the amplitude of u, thereby making the scaling of u inappropriate. 


In case I = V = 0 so that a = 6 =0, (2.72) predicts (A= m = 1, w= 87}, 
y=1) 


59 
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u(#) = (6-2 -1)~2sin (50 = =) sin (50+) , 


which has an amplitude about 2 for 6 >> 1. Figure 2.10 shows a case. 


dt=0.0392699 


Fig. 2.10 Forced undamped vibrations with rapid forcing. 


With a = 0.058? = 5, we get a significant contribution from the free vibra- 
tions (the homogeneous solution of the ODE) as shown in Figure 2.11. For 
larger a values, one must base u, on J instead. (The graphs in Figure 2.10 
and 2.11 were produced by numerical simulations with 160 time steps per 
period of the forcing.) 


Displacement scale based on J. Choosing uc = I gives 


——+u=yc0s(6t), u(0)=1, w(0)=8, (2.79) 


with 
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dt=0.0392699 


Fig. 2.11 Forced undamped vibrations with rapid forcing and initial displacement of 
5. 


Vte V jm 

e NT a 
t24 A A 

pe Se (2.81) 


This scaling is not relevant close to resonance since then uc > I. 


2.2.4 Damped vibrations with forcing 


We now introduce a linear damping force bu’(t) in the equation of motion: 


mu” +bu' + ku = Acos(vt), u(0)=J, u/(0)=V. (2.82) 


Figure 2.12 shows a typical one-degree-of-freedom mechanical system with a 
linear dashpot, representing the damper (bu’), a linear spring (ku), and an 
external force (F). 

The standard scaling procedure results in 
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b>) 


Fig. 2.12 Oscillating body with external force, attached to a spring and damper. 


üu tebdūù tk t2 I 
B: a te u= —— Acos(Ytt), u(0)=—, u (0)= ri : 
dt? mdt m Muc Uc Uc 


(2.83) 


The exact solution. As always, it is a great advantage to look into exact 
solutions for controlling our choice of scales. Using SymPy to solve (2.82) is, 
in principle, very straightforward: 


>>> diffeq = diff(u(t), t, t) + b/m*diff (u(t), t) + we*2*u(t) 
>>> s = dsolve(diffeq, u(t)) 

>>> s.rhs 

Cixexp(t*(-b - sqrt(b - 2*m*w)*sqrt(b + 2*m*w))/(2*m)) + 
C2*exp(t*(-b + sqrt(b - 2*m*w)*sqrt(b + 2*m*w))/(2*m)) 


This is indeed the correct solution, but it is on a complex exponential func- 
tion form, valid for all b, m, and w. We are interested in the case with 
small damping, b < 2mw, where the solution is an exponentially damped 
sinusoidal function. Rewriting the expression in the right form is tricky with 
SymPy commands. Instead, we demonstrate a common technique when doing 
symbolic computing: general procedures like dsolve are replaced by manual 
steps. That is, we solve the ODE “by hand”, but use SymPy to assist the 
calculations. 

The solution is composed of a homogeneous solution up of mu” + bu! + 
ku = 0 and one particular solution up of the nonhomogeneous equation 
mu” + bu! + ku = Acos(wt). The homogeneous solution with damped oscil- 
lations (requiring b < 2V/mk) can be found by the following code. We have 
divided the differential equation by m and introduced B = 5b/ m and let A1 
represent A/m to simplify expressions and help SymPy with less symbols 
in the equation. Without these simplifications, SymPy stalls in the compu- 
tations due to too many symbols in the equation. The problem is actually 
a solid argument for scaling differential equations before asking SymPy to 
solve them since scaling effectively reduces the number of parameters in the 
equations! 

The following SymPy steps derives the solution of the homogeneous ODE: 
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u = symbols(’u’, cls=Function) 
t, w, B, A, Al, m, psi = symbols(’t w B A Al m psi’, 
positive=True, real=True) 


def ode(u, homogeneous=True) : 
h = diff(u, t, t) + 2*Bediff(u, t) + w**2*u 
f = Alxcos(psi*t) 
return h if homogeneous else h - f 


# Find coefficients in polynomial (in r) for exp(r*t) ansatz 
r = symbols(’r’) 

ansatz = exp(r*t) 

poly = simplify (ode(ansatz) /ansatz) 


# Convert to polynomial to extract coefficients 
poly = Poly(poly, r) 

# Extract coefficients in poly: a_*t**2 + b_*t + c_ 
a_, b_, c_ = poly.coeffs() 

# Assume b_**2 - 4*a_*c_ < 0 

d = -b_/(2*a_) 

ifa ==" Ie: 


omega = sqrt(c_ - (b_/2)**2) # nicer formula 
else: 
omega = sqrt (4*a_*c_ - b_**2)/(2*a_) 


# The homogeneous solution is a linear combination of a 
# cos term (u1) and a sin term (u2) 

ul = exp(d*t) *cos(omega*t) 

u2 = exp(d*t) *sin(omega*t) 

C1, C2, V, I = symbols(’C1 C2 V I’, real=True) 

u_h = simplify(Ci*u1 + C2*u2) 

print evens uh 


The print out shows 


un = e7 Pt (c cos( y w? — B?t) + C2sin( V w? — B?t)) ; 


where C1 and C2 must be determined by the initial conditions later. It is 
wise to check that up is indeed a solution of the homogeneous differential 
equation: 


assert simplify(ode(u_h)) == 


We have previously just printed the residuals of the ODE and initial condi- 
tions after inserting the solution, but it is better in a code to let the program- 
ming language test that the residuals are symbolically zero. This is achieved 
using the assert statement in Python. The argument is a boolean expres- 
sion, and if the expression evaluates to False, an AssertionError is raised 
and the program aborts (otherwise assert runs silently for a True boolean 
expression). Hereafter, we will use assert for consistency checks in computer 
code. 
The ansatz for the particular solution up is 


Up = C3 cos(Yt) + Cysin(yt), 
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which inserted in the ODE gives two equations for C3 and C4. The relevant 
SymPy statements are 


# Particular solution 

C3, C4 = symbols(’C3 C4’) 

u_p = C3*cos(psitt) + C4*sin(psi*t) 

eqs = simplify(ode(u_p, homogeneous=False) ) 


# Collect cos(omega*t) terms 

print ’eqs:’, eqs 

eq_cos = simplify(eqs.subs(sin(psi*t), 0).subs(cos(psi*t), 1)) 
eq_sin = simplify(eqs.subs(cos(psi*t), 0).subs(sin(psi*t), 1)) 
s = solve([eq_cos, eq_sin], [C3, C4]) 

u_p = simplify(u_p.subs(C3, s[C3]).subs(C4, s[C4])) 


# Check that the solution is correct 
assert simplify(ode(u_p, homogeneous=False)) == 


Using the initial conditions for the complete solution u = up + up determines 
Cı and Co: 


u_sol = uh + u_p # total solution 

# Initial conditions 

eqs = [u_sol.subs(t, 0) - I, u_sol.diff(t).subs(t, 0) - V] 
# Determine Ci and C2 from the initial conditions 

s = solve(eqs, [C1, C2]) 

u_sol = u_sol.subs(Ci, s[Ci]).subs(C2, s[C2]) 


Finally, we should check that u_sol is indeed the correct solution: 


checks = dict( 
ODE=simplify(expand(ode(u_sol, homogeneous=False))), 
IC1=simplify(u_sol.subs(t, 0) - I), 
IC2=simplify(diff(u_sol, t).subs(t, 0) - V)) 

for check in checks: 
msg = ’/s residual: %s’ % (check, checks[check] ) 
assert checks[check] == sympify(0), msg 


Finally, we may take u_sol = u_sol.subs(A, A/m) to get the right expres- 
sion for the solution. Using latex (u_sol) results in a huge expression, which 
should be manually ordered to something like the following: 


2 Am7! 
— 4B? + 2? 


= (cy00s(¢ w? — B?) + C2sin (t 2B) ) 


Am-!Q4 41 By? + 1:9? 
~ 4B2yy2 + N2 
Cy = — Amt BR +4I B3y? + IBN? + 4V BY? +. Vin? 
Vw? — B? (4B2y?2 + 2?) 


u (2Bysin (wt) — R cos (yt)) + 


Cı 


Q=Y-w?. 
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The most important feature of this solution is that there are two time 
scales with frequencies 7 and Vw? — B?, respectively, but the latter appears 
in terms that decay as e~ P+ in time. The attention is usually on longer periods 
of time, so in that case the solution simplifies to 


A —1 
i= e (2Bysin (Wt) — 2 cos (Wt) 
A 1 n (pw) t 
= in an r O yaya 


= for (1 +Q?(8-871)) 7? cos(Yt+ ¢), (2.84) 


where we have introduced the dimensionless numbers 


and 


= 2B = Q7! 
ġ = tan 1 (=) = tan 1 (2) š 


Q is commonly called quality factor and ¢ is the phase shift. Dividing (2.84) 
by A/k, which is a common scale for u, gives the dimensionless relation 


FE = SERQA) costut +0), R(Q,ô) = (14+.Q2(6-671))*. (2.85) 


Choosing scales. Much of the discussion about scales in the previous sec- 
tions are relevant also when damping is included. Although the oscillations 
with frequency Vw? — B2 die out for t> BT}, we start with using this fre- 
quency for the time scale. A highly relevant assumption for engineering ap- 
plications of (2.82) is that the damping is small. Therefore, Vw? — B? is close 
to w and we simply apply te = 1/w as before (if not the interest in large t for 
which the oscillations with frequency w has died out). 

The coefficient in front of the u’ term then becomes 


The rest of the ODE is given in the previous section, and the particular 
formulas depend on the choices of te and Ue. 


Choice of u, at resonance. The relevant scale for u, at or nearby resonance 
(Y = w) becomes different from the previous section, since with damping, the 
maximum amplitude is a finite value. For t> B71, when the sinwt term is 
dominating, we have for Y = w: 
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Am—!2Bw A 
~ 4B? sin(wt) = 2Bm 2Bmd sin(wt) = by a ees 
This motivates the choice 
7 A _ A 
te ib bw 


(It is wise during computations like this to stop and check the dimensions: 
A must be [MLT~?] from the original equation (F(t) must have the same 
dimension as mu”), bu’ must also have dimension [MLT~?], implying that 
b has dimension [MT~*]. A/b then has dimension LTT}, and A/(bw) gets 
dimension [|L], which matches what we want for uc.) 

The differential equation on dimensionless form becomes 


G i — 7 tus = ycos(ôt), u(0)=a, w(0)=8, (2.86) 
with 
I Ib Jk 
Vt. Vb 
B= ae (2.88) 
2 
ead (2.89) 
Muc k 
te 4% 
ô T (2.90) 


Choice of u, when w > w. In the limit w > w and t> B7}, 
A A 
u x —; cospt = — cos yt, 
mw k 


showing that uc = A/k is an appropriate displacement scale. (Alternatively, 
we get this scale also from demanding y = 1 in the ODE.) The dimensionless 
numbers a, 3, and 6 are as for the forced vibrations without damping. 


Choice of uce when w & y. In the limit w < yY, we should base te on the 
rapid variations in the excitation: te = 1/4. 


Software. It is easy to reuse a solver for a general vibration problem also 
in the dimensionless case. In particular, we may use the solver function in 
the file vib. py: 


def solver(I, V, m, b, s, F, dt, T, damping=’linear’): 


for solving the ODE problem 
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mu” + f(u')+s(u) = F(t), u(0) =I, u'(0) =V, te (0,7), 


with time steps dt. With damping=’ linear’, we have f(u’) = bu’, while the 
other value is ’quadratic’, meaning f(u’) = b|u’|u’. Given the dimensionless 
numbers a, 6, y, 6, and Q, an appropriate call for solving (2.73) is 


u, t = solver(I=alpha, V=beta, m=1, b=1.0/Q, 
s=lambda u: u, F=lambda t: gamma*cos(delta*t), 
dt=2*pi/n, T=2*pixP) 


where n is the number of intervals per period and P is the number of periods 
to be simulated. We way wrap this call in a solver_scaled function and 
wrap it furthermore with joblib to avoid repeated calls, as we explained in 
Section 2.1.4: 


from vib import solver as solver_unscaled 


def solver_scaled(alpha, beta, gamma, delta, Q, T, dt): 


Solve u’? + (1/Q)*u’ + u = gamma*cos(delta*t) , 
u(0)=alpha, u’(1)=beta, for (0,T] with step dt. 


print ’Computing the numerical solution’ 

from math import cos 

return solver_unscaled(I=alpha, V=beta, m=1, b=1./Q, 
s=lambda u: u, 
F=lambda t: gamma*cos(delta*t), 
dt=dt, T=T, damping=’ linear’) 


import joblib 
disk_memory = joblib.Memory(cachedir=’ temp’) 
solver_scaled = disk_memory.cache(solver_scaled) 


This code is found in vib_scaled.py and features an application for running 
the scaled problem with options on the command-line for a, 6, y, 6, Q, num- 
ber of time steps per period, and number of periods (see the main function). 
It is an ideal application for exploring scaled vibration models. 


2.2.5 Oscillating electric circuits 


The differential equation for an oscillating electric circuit is very similar to the 
equation for forced, damped, mechanical vibrations, and their dimensionless 
form is identical. This fact will now be demonstrated. 

The current I(t) in a circuit having an inductor with inductance L, a 
capacitor with capacitance C, and overall resistance R, obeys the equation 


a Re 1 
5 lee = VO: (2.91) 


where V(t) is the voltage source powering the circuit. We introduce 
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and get 
PI tRd ie Ja 
d L d EO k 
Here, we have scaled V(t) according to 
i V (tet) 
V(t) = —— =. 
() max; V (t) 


The time scale te is chosen to make Ï and I/(LC) balance, te = VLC. 
Choosing Ie to make the coefficient in the source term of unit size, means 
I, = LCV... With 


2 — 
teVe V(t). 


C 
-1_ La 
Q= R7, 
we get the scaled equation 
dI yar Sas 
— ~ —+I=V(t 2.92 
tE +I =V, (2.92) 


which is basically the same as we derived for mechanical vibrations. (Two 
additional dimensionless variables will arise from the initial conditions for J, 
just as in the mechanics cases.) 
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Chapter 3 


Basic partial differential equation 
models 


This chapter extends the scaling technique to well-known partial differential 
equation (PDE) models for waves, diffusion, and transport. We start out with 
the simplest 1D models of the PDEs and then progress with additional terms, 
different types of boundary and initial conditions, and generalizations to 2D 
and 3D. 


3.1 The wave equation 


A standard, linear, one-dimensional wave equation problem in a homogeneous 
medium may be written as 


Pu Au 
I T Bax? 
where c is the constant wave velocity of the medium. With a briefer notation, 
where subscripts indicate derivatives, the PDE (3.1) can be written ut = 
Cuzz. This subscript notation will occasionally be used later. 
For any number of dimensions in heterogeneous media we have the gener- 
alization 


x €(0,L), te (0,7), (3.1) 


2 
CH HaV (PVu)+f, my2€0, te (0T), (3.2) 


where f represents a forcing. 


3.1.1 Homogeneous Dirichlet conditions in 1D 


Let us first start with (3.1), homogeneous Dirichlet conditions in space, and 
no initial velocity uz: 
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u(x,0) = I(x), x € [0, L], (3.3) 
o 
gee =0, x € (0, L], (3.4) 
u(0,t) =0, te (0,7), (3.5) 
u(L,t) =0, te (0,7). (3.6) 


The independent variables are x and t, while u is the dependent variable. 
The rest of the parameters, c, L, T, and I(x), are given data. 

We start with introducing dimensionless versions of the independent and 
dependent variables: 


C= yp GS UES . 
Le te Uc 


Inserting the z = zeg, etc., in (3.1) and (3.3)-(3.6) gives 


Pu Peou 


P2 = g2 ðr? TE (0,L/ze), te (0,T/tel, 
a(z, 0) = EP., Z€ [0,L/zel, 
Uc 
O_ E 
gE =0, zE [0,L/z£e], 
u(0,t) =0, te (0,T/tel, 
u(L/ze,t) =0, te (0,T/te. 


The key question is how to define the scales. A natural choice is £e = L 
since this makes Z € [0,1]. For the spatial scale and the problem governed by 
(3.1) we have some analytical insight that can help. The solution behaves like 


u(a,t) = fr(w—ct)+ fr(at+ct), (3.7) 


i.e., a right- and left-going wave with velocity c. The initial conditions con- 
strain the choices of fr and fr to ft + fr =I and —cf;, +cfp =0. The 
solution is fR = fr = 5 and consequently 


1 1 
u(x,t) = gi le — et) + z1 0+ et), 


which tells that the initial condition splits in two, half of it moves to the 
left and half to the right. This means in particular that we can choose 
Ue = Max, |I (x)| and get |u| < 1, which is a goal. It must be added that 
boundary conditions may result in reflected waves, and the solution is then 
more complicated than indicated in the formula above. 
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Regarding the time scale, we may look at the two terms in the scaled PDE 
and argue that if |u| and its derivatives are to be of order unity, then the size 
of the second-order derivatives should be the same, and te can be chosen to 
make the coefficient t2c?/x? unity, i.e., te = L/c. Another reasoning may set 
te as the time it takes the wave to travel through the domain [0, L]. Since the 
wave has constant speed c, te = L/c. 

With the described choices of scales, we end up with the dimensionless 
initial-boundary value problem 


aes ze (0,1), te (0,T], (3.8) 
ue) oan z€ (0,1), (3.9) 
Ža(z.0) =0, z € [0,1], (3.10) 
u(0,£) =0, te (0,7), (3.11) 
u(1,t) =0, te (0,7). (3.12) 


Here, T = Tc/L. 

The striking feature of (3.8)-(3.12) is that there are no physical parameters 
involved! Everything we need to specify is the shape of the initial condition 
and then scale it such that it is less than or equal to 1. 

The physical solution with dimension is recovered from w(Z,t) through 


u(x,t) = oe I(x) u(ZL,tL/c) (3.13) 


3.1.2 Implementation of the scaled wave equation 


How do we implement (3.8)-(3.12)? As for the simpler mathematical models, 
we suggest to implement the model with dimensions and observe how to set 
parameters to obtain the scaled model. In the present case, one must choose 
L=1, c=1, and scale I by its maximum value. That’s all! 

Several implementations of 1D wave equation models with different degree 
of mathematical and software complexity come along with these notes. The 
simplest version is wave1D_u0.py that implements (3.1) and (3.3)-(3.6). This 
is the code to be used in the following. It is described in Section 2.3 in [7]. 
Waves on a string. As an example, we may let the original initial-boundary 
value problem (3.1)-(3.6) model vibrations of a string on a string instrument 
(e.g., a guitar). With u as the displacement of the string, the boundary con- 
ditions u = 0 at the ends are relevant, as well as the zero velocity condition 
Ou/Ot =0 at t=0. The initial condition I(x) typically has a triangular shape 
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for a picked guitar string. The physical problem needs parameters for the am- 
plitude of I(x), the length L of the string, and the value of c for the string. 
Only the latter is challenging as it involves relating c to the pitch (i.e., time 
frequency) of the string. In the scaled problem, we can forget about all this. 
We simply set L= 1, c= 1, and let I(x) have a peak of unity at x = xo € (0,1): 


I(x) a ie £ < T0, 
max, I(x) (1—2)/(1—20), otherwise 


The dimensionless coordinate of the peak, xo, is the only dimensionless pa- 
rameter in the problem. For fixed xo, one single simulation will capture all 
possible solutions with such an initial triangular shape. 


Detecting an already computed case. The file wave1D_u0_scaled.py 
has functionality for detecting whether a simulation corresponds to a previ- 
ously run scaled case, and if so, the solution is retrieved from file. The im- 
plementation technique makes use of joblib, but is more complicated than 
shown previously in these notes since some of the arguments to the function 
that computes the solution are functions, and one must recognized if the 
function has been used as argument before or not. There is documentation 
in the wave1D_u0_scaled.py file explaining how this is done. 


3.1.3 Time-dependent Dirichlet condition 


A generalization of (3.1)-(3.6) is to allow for a time-dependent Dirichlet con- 
dition at one end, say u(0,t) = Uz(t). At the other end we may still have 
u=0. This new condition at x = 0 may model a specified wave that enters 
the domain. For example, if we feed in a monochromatic wave Asin(k(a—ct)) 
from the left end, Uz (t) = Asin(kct). This forcing of the wave motion has its 
own amplitude and time scale that could affect the choice of ue and te. 

The main difference from the previous initial-boundary value problem is 
the condition at x = 0, which now reads 


u(0,t) = Ur (tte) 


Uc 


in scaled form. 


Scaling. Regarding the characteristic time scale, it is natural to base this 
scale on the wave propagation velocity, together with the length scale, and not 
on the time scale of Uz (t), because the time scale of Uz basically determines 
whether short or long waves are fed in at the boundary. All waves, long 
or short, propagate with the same velocity c. We therefore continue to use 
te = L/c. 

The solution u will have one wave contribution from the initial condition I 
and one from the feeding of waves at x = 0. This gives us three choices of ue: 
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max, |I| + max;|Uz|, max, |I|, or max: |Uz|. The first seems relevant if the 
size of I and Uy are about the same, but then we can choose either max, |T| 
or max; |Uz| as characteristic size of u since a factor of 2 is not important. If 
I is much less than Uz, uc = max; |uz| is relevant, while uc = max, |I| is the 
choice when J has much bigger impact than Uz on u. 

With ue = max; |Uz(t)|, we get the scaled problem 


eee ze (0,1), te (0,7), (3.14) 
u(#,0) = ay ze [0,1], (3.15) 

© (2,0) = 0, z € [0,1], (3.16) 
u(0,t) = men te (0,7], (3.17) 
u(1,t) =0, te (0,7). (3.18) 


Also this problem is free of physical parameters like c and L. The input is 
completely specified by the shape of I(x) and Uz/(t). 

Software. Software for the original problem with dimensions can be reused 
for (3.14)-(3.18) by setting L = 1, c=1, and scaling Uz(t) and I(x) by 
max; |Uz(t)|. 


Specific case. As an example, consider 


U,(t) =asin(wt) forO<t< 2, else 0, 
T 
Ha) = Ae- (@-1/2?/0 , 


That is, we start with a Gaussian peak-shaped wave in the center of the 
domain and feed in a sinusoidal wave at the left end for two periods. The 
solution will be the sum of three waves: two parts from the initial condition, 
plus the wave fed in from the left. 

Since max; |UL| =a we get 


ii(Z,0) = te E/2)? 3)? (3.19) 
a 
u(0,t) =sin(twL/c). (3.20) 


Here, Uz, models an incoming wave asin(k(x —ct), with k specified. The re- 
sult is incoming waves of length A = 27/k. Since w = kc, u(0,t) = sin(kLt) = 
sin(27tL/X). (This formula demonstrates the previous assertion that the time 
scale of Uz, i.e., 1/w, determines the wave length 1/w = A/(27) in space.) We 
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realize from the formulas (3.19) and (3.20) that there are three key dimen- 
sionless parameters related to these specific choices of initial and boundary 
conditions: 
A L L 
Q = —, B=-, y=kL=2n-. 
a o À 


With a, 8, and y we can write the dimensionless initial and boundary con- 
ditions as 


275 2 
u(Z,0) = ae? (@-3) x 
u(0,t) = sin(yt). 
The dimensionless parameters have the following interpretations: 


e a: ratio of initial condition amplitude and amplitude of incoming wave at 
xr=0 

e £: ratio of length of domain and width of initial condition 

e ~q: ratio of length of domain and wave length of incoming wave 


Again, these dimensionless parameters tell a lot about the interplay of the 
physical effects in the problem. And only some ratios count! 
We can simulate two special cases: 


1. a= 10 (large) where the incoming wave is small and the solution is dom- 
inated by the two waves arising from I(x), 

2. a=0.1 (small) where the incoming waves dominate and the solution has 
the initial condition just as a small perturbation of the wave shape. 


We may choose a peak-shaped initial condition: 8 = 10, and also a relatively 
short incoming wave compared to the domain size: y = 6r (i.e., wave length of 
incoming wave is L/6). A function simulate_Gaussian_and_incoming_ wave 
in the file session.py applies the general unscaled solver in wave1D_dn. 
py for solving the wave equation with constant c, and any time-dependent 
function or 0u/Ox = 0 at the end points. This solver is trivially adapted to 
the present case. Figures 3.1 and 3.2 shows snapshots of how u(Z,t) evolves 
due to a large/small initial condition and small/large incoming wave at the 
left boundary. 


Movie 1: a= 10. https://github.com/hplgit/scaling-book/raw/master/doc/ 
pub/book/htm1/mov-scaling/gaussian_plus_incoming/alpha10.mp4 


Movie 2: a=0.1. https: //github.com/hplgit/scaling-book/raw/master/doc/ 
pub/book/htm1/mov-scaling/gaussian_plus_incoming/alpha01.mp4 
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t=0.000, u(0,t)=U_O(t), du(L,t)/dx=0 t=0.250, u(0,t)=U_O(t), du(L,t)/dx=0 
10 10 
5 s 
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Fig. 3.1 Snapshots of solution with large initial condition and small incoming wave 
(a = 10). 


3.1.4 Velocity initial condition 


Now we change the initial condition from u = I and 0u/0t = 0 to 


u(x,0)=0, (3.21) 


a) 
gE =V (g). (3.22) 


Impact problems are often of this kind. The scaled version of u;(2,0) = V(x) 
becomes 


Analytical insight. From (3.7) we now get fz + fr =0 and cf; —cfp=V. 
Introducing W(x) such that W'(x) = V(x), a solution is fg = 5W/c and 
—fr= 5W/ c. We can express this solution through the formula 
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t=0.000, u(0,t)=U_O(t), du(L,t)/dx=0 _t=0.250, u(0,t)=U_O(t), du(L,t)/dx=0 
1.0 1.0 
05 05 l 
bolion T 3.68 | — 
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x x 
t=0.500, u(0,t)=U_O(t), du(L,t)/dx=0 t=1.375, u(0,t)=U_O(t), du(L,t)/dx=0 
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Fig. 3.2 Snapshots of solution with small initial condition and large incoming wave 
(a=0.1). 


1 xz+ct 1 
E f ve- Weta)—Wes ai. (3.23) 
20 Jaca 2c 
Scaling. Since V is the time-derivative of u, the characteristic size of V, call 
it Vo, is typically uc/te. If we, as usual, base te on the wave speed, te = L/c, 
we get uc = V.L/c. Looking at the solution (3.23), we see that ue has size 
mean(V)L£/(2c), where mean(V) is the mean value of V (W ~ mean(V)L). 
This result suggests V. = mean(V) and ue = mean(V)L/(2c). One may argue 
that the factor 2 is not important, but if we want |u| € [0,1] it is convenient 
to keep it. 
The scaled initial condition becomes 


Nonzero initial shape. Suppose we change the initial condition u(x,0) =0 
to u(x,0) = I(x). The scaled version of this condition with the above ue based 
on V becomes 


2cI (Txe) 


a(z,0) = Lmean(V) (3.24) 
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Check that dimensionless numbers are dimensionless! 


Is a dimensionless number really dimensionless? It is easy to make er- 
rors when scaling equations, so checking that such fractions are dimen- 
sionless is wise. The dimension of I is the same as u, here taken to be 
displacement: [L]. Since V is ðu/ðt, its dimension is [LT~*]. The dimen- 
sions of c and L are [LT~'] and [L]. The dimension of the right-hand 
side of (3.24) is then 


[LTH _ 
pe 


demonstrating that the fraction is indeed dimensionless. 


F 


One may introduce a dimensionless initial shape, T(z) = I(@L)/ maxa |I]. 
Then 


u(z,0) = al (z), 
where a the dimensionless number 


_ 2cmax,|I(2)| 
= L mean(V) 


If V is much larger than J, one expects that the influence of I is small. 
However, it takes time for the initial velocity V to influence the wave motion, 
so the speed of the waves c and the length of the domain L also play a 
role. This is reflected in a, which is the important parameter. Again, the 
scaling and the resulting dimensionless parameter(s) teach us much about 
the interaction of the various physical effects. 


3.1.5 Variable wave velocity and forcing 


The next generalization regards wave propagation in a non-homogeneous 
medium where the wave velocity c depends on the spatial position: c= c(x). 
To simplify the notation we introduce A(x) = c?(x). We introduce homoge- 
neous Neumann conditions at x =0 and x = L. In addition, we add a force 
term f(x,t) to the PDE, modeling wave generation in the interior of the 
domain. For example, a moving slide at the bottom of a fjord will generate 
surface waves and is modeled by such an f(x,t) term (provided the length of 
the waves is much larger than the depth so that a simple wave equation like 
(3.25) applies). The initial-boundary value problem can be then expressed as 


78 3 Basic partial differential equation models 


u ð Ou 
= t L), t T 2 
TE AWZ) — ee OL), t€O.T}, 6.25) 
u(z,0) = I(x), xe [0, L], (3.26) 
© (2,0) =0, xe f0, L], (3.27) 
208) =0, te(0,T], (3.28) 
Ox 
LD =0, te(0,T]. (3.29) 
Ox 
Non-dimensionalization. We make the coefficient A non-dimensional by 
\) = Í (3.30) 


where one normally chooses the characteristic size of A, Ac, to be the maxi- 
mum value such that |A| < 1: 


Ac= max A(x). 
x€(0,L) 


Similarly, f has a scaled version 


f. F(Exc,tte) 


f(z, = r 


where normally we choose 
fe = max|f(x,t)|. 
x,t 


Inserting dependent and independent variables expressed by their non- 
dimensional counterparts yields 


Pu HO (> en tote = a = 2 
or aes GEZ) a, e01), Fe (0.7, 
cos le z € [0,1], 
Uc 
IR =0, z € [0,1], 
2 50,8 =0, te (0,T], 
Ca E te (0,T], 
Ox 


with T = Tc/L. 
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Choosing the time scale. The time scale is, as before, chosen as te = 
L/V\e. Note that the previous (constant) wave velocity c now corresponds 
to \/X(x). Therefore, V/A; is a characteristic wave velocity. 

One could wonder if the time scale of the force term, f(x,t), should influ- 
ence te, but as we reasoned for the boundary condition u(0,t) = Uz (t), we let 
the characteristic time be governed by the signal speed in the medium, i.e., 
by VX¢ here and not by the time scale of the excitation f, which dictates the 
length of the generated waves and not their propagation speed. 


Choosing the spatial scale. We may choose ue as max, |I(x)|, as before, 
or we may fit we such that the coefficient in the source term is unity, i.e., all 
terms balance each other. This latter idea leads to 


L? fe 
Urp = 
c A: 


and a PDE without parameters, 


Pu əf- əN =_- 
aE = 5 Ow e) + FED. 


T 
The initial condition u(x,0) = I(x) becomes in dimensionless form 


where 


In the case ue = max, |I (x)|, U(z,0) = I(Z) and the 6 parameter appears 
in the PDE instead: 


25 A u D 

Pa 2 (un Bt) eaaa 
With V = 0, and u =0 or uz = 0 on the boundaries x = 0, L, this scaling 
normally gives |u| < 1, since initially |Z| < 1, and no boundary condition can 
increase the amplitude. However, the forcing, f, may inherit spatial and tem- 
poral scales of its own that may complicate the matter. The forcing may, for 
instance, be some disturbance moving with a velocity close to the propaga- 
tion velocity of the free waves. This will have an effect akin to the resonance 
for the vibration problem discussed in section 2.2.2 and the waves produced 
by the forcing may be much larger than indicated by 8. On the other hand, 
the forcing may also consist of alternating positive and negative parts (ret- 
rogressive slides constitute an example). These may interfere to reduce the 
wave generation by an order of magnitude. 


Scaling the velocity initial condition. The initial condition u;(x,0) = 
V(x) has its dimensionless variant as 
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=. te V(Lī) 
NS ee WO 
which becomes 
OU ,_ L maxs |V (x)| z ner. o 
py 9) = x aa V(z), if Ue = max |I (x)|, 
or 
= 2 
M gja VE mO a a aa ymax fet) 


Ot L maxz+|f(x,t)| 


Introducing the dimensionless number a (cf. Section 3.1.4), 


gat We_maxe|V(2)| 


© L maxg t| f(c, e| 


we can write 


— 

R 
© 

© 
II 


OU ,_ a~1V(z), Ue = Max, |T| 
aa BTV (3), Uc = te fe 
3.1.6 Damped wave equation 


A linear damping term b0u/Ot is often added to the wave equation to model 
energy dissipation and amplitude reduction. Our PDE then reads 


3u ðu 0 Ou 
E toz a7 O- ) +e t). (3.31) 


The scaled equation becomes 


ü tedn tds 3 (~, 3T fez n 
db Of L? OF Owe) Pa 


The damping term is usually much smaller than the two other terms involv- 
ing u. The time scale is therefore chosen as in the undamped case, te = L/ VAc. 
As in Section 3.1.5, we have two choices of ue: Ue = max, |I| or Uc = t2 fe. 
The former choice of ue gives a PDE with two dimensionless numbers, 


Uc 


Pu ðu ð Ou 
158 = i AOS) +879, (3.32) 


where 
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Y ST 


measures the size of the damping, and £ is as given in Section 3.1.5. With 
uc = t? fe we get a PDE where only y enters, 


u ðu ðf ôu) + 
= Tt == A(T + f (x,t). 3.33 
aE tS = be OW ge) +F@D (3.33) 
The scaled initial conditions are as in Section 3.1.5, so in this latter case 8 
appears in the initial condition for u. 
To summarize, the effects of V, f, and damping are reflected in the di- 
mensionless numbers a, 3, and y, respectively. 


3.1.7 A three-dimensional wave equation problem 


To demonstrate how the scaling extends to in three spatial dimensions, we 


consider 
Oru 0 Ou 0 Ou o Ou 
— = ; .34 
OP” Oa OZ) OZ) 0%) (339 
Introducing 
_ x — y a 22 af Z u 
t= —, y= —, Z= —, t=—, u= —, 
Te Ye Ze te Ue 


and scaling À as À = A(Z£e, 9 Yc, ZZ) /Ae, we get 
Ot _ tere ə (59%) | tere ð (5%). Ac ð (525). 
Ot? x2 02 \ OF y2 Oy \ OY z2 ðz \ OZ 
Often, we will set £e = ye = Ze = L where L is some characteristic size of the 


domain. As before, te = L/VAc, and these choices lead to a dimensionless 
wave equation without physical parameters: 


Ou A (du ð (-0u ð (-0u 
OP Oz OZ) ag (az) Ta OZ) oe) 


The initial conditions remain the same as in the previous one-dimensional 
examples. 


3.2 The diffusion equation 


The diffusion equation in a one-dimensional homogeneous medium reads 
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ðu _ Ou 

at Ox?’ 
where a is the diffusion coefficient. The multi-dimensional generalization to 
a heterogeneous medium and a source term takes the form 


x € (0, L), t€ (0,7), (3.36) 


Ĉu L Y.(aVu) +f, my2€2, te(0T]: (3.37) 


We first look at scaling of the PDE itself, and thereafter we discuss some 
types of boundary conditions and how to scale the complete initial-boundary 
value problem. 


3.2.1 Homogeneous 1D diffusion equation 


Choosing the time scale. To make (3.36) dimensionless, we introduce, as 
usual, dimensionless dependent and independent variables: 
x - ¢ = u 
; t=) “U= ; 
Le te Ue 


Inserting the dimensionless quantities in the one-dimensional PDE (3.36) 
results in 


du tea Pu 7 x 
2S T 1), t T=T/te]. 
Ot L2 0x2’ x (0, ), € (0, / c] 
Arguing, as for the wave equation, that the scaling should result in 
Ou ru 
ar oR 


of the same size (about unity), implies t.a/L? = 1 and therefore te = L? /a. 


Analytical insight. The best way to obtain the scales inherent in a problem 
is to obtain an exact analytic solution, as we have done in many of the ODE 
examples in this booklet. However, as a rule this is not possible. Still, often 
highly simplified analytic solutions can be found for parts of the problem, or 
for some closely related problem. Such solutions may provide crucial guidance 
to the nature of the complete solution and to the appropriate scaling of the 
full problem. We will employ such solutions now to learn about scales in 
diffusion problems. 

One can show that u = Ae~?‘sin(kz) is a solution of (3.36) if p= ak?, for 
any k. This is the typical solution arising from separation of variables and 
reflects the dynamics of the space and time in the PDE. Exponential decay 
in time is a characteristic feature of diffusion processes, and the e-folding 
time can then be taken as a time scale. This means te = 1/p ~ k~?. Since k 
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is related to the spatial wave length A through k = 27/X, it means that te 
depends strongly on the wave length of the sine term sin(ka). In particular, 
short waves (as found in noisy signals) with large k decay very rapidly. For 
the overall solution we are interested in how the longest meaningful wave 
decays and use that time scale for te. The longest wave typically has half a 
wave length over the domain [0, L]: u = Ae~” sin(z2/L) (k =7/L), provided 
u(0,t) = u(L,t) =0 (with uz(L,t) =0, the longest wave is L/4, but we look 
at the case with the wave length L/2). Then te = L?/an~?, but the factor 
mT? is not important and we simply choose te = L?/a, which equals the time 
scale we arrived at above. We may say that te is the time it takes for the 
diffusion to significantly change the solution in the entire domain. 

Another fundamental solution of the diffusion equation is the diffusion of a 
Gaussian function: u(x,t) = K(4mat)~!/? exp (—2x?/(4at)), for some constant 
K with the same dimension as u. For the diffusion to be significant at a 
distance x = L, we may demand the exponential factor to have a value of 
e71 ~ 0.37, which implies t = L?/(4a), but the factor 4 is not of importance, 
so again, a relevant time scale is te = L?/a. 


Choosing other scales. The scale ue is chosen according to the initial 
condition: ue = max,¢0,z) | (x)|. For a diffusion equation uz = aug, with u = 
0 at the boundaries x = 0, L, the solution is bounded by the initial condition 
I(x). Therefore, the listed choice of ue implies that |u| < 1. (The solution 
u = Ae~”'sin(kx) is such an example if k = nr/L for integer n such that 
u=0 for z = 0 and z = L.) 

The resulting dimensionless PDE becomes 


ðu = Ou z = 
— = oe 1), t T i 
a = pg? FE G1), te (0,7), (3.38) 
with initial condition 
ae sai I(x) 
= [ — ere 


Notice that (3.38) is without physical parameters, but there may be param- 
eters in I(x). 


3.2.2 Generalized diffusion PDE 


Turning the attention to (3.37), we introduce the dimensionless diffusion 
coefficient 


a(z,¥,Z) = az talte, Yo, Zc2), 


typically with 
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Qe = maxa(x,y,Z). 
XLyY,z 


The length scales are 


We scale f in a similar fashion: 


(E921) = fo" S(Exe,Gyc%e, tte), 
with 
fe = max \f(x,y,z,t)| i 
LY, zt 
Also assuming that £e = yo = Ze = L, and ue = maxz,y (I(x, y,z), we end up 


with the scaled PDE 


Ou a _ = _ 

Sp =Y (ava) + Bf, 2,9,2¢ Q, te (0,7). (3.39) 
Here, V means differentiation with respect to dimensionless coordinates 7, 
y, and Z. The dimensionless parameter (3 takes the form 


fe tefe _ L? maXr y,z,t |f(x,y,z,t)| 
Uc Qa MaxXy yz |I (x,y, z)| 


B 


The scaled initial condition is ū = Í as in the 1D case. 
An alternative choice of uc is to make the coefficient t.f-/uce in the source 
term unity. The scaled PDE now becomes 
due 2 
=Y (ava) +i, (3.40) 
but the initial condition features the 6 parameter: 


The 6 parameter can be interpreted as the ratio of the source term and 
the terms with u: 


kaU ty cs oy os lel 
uc/te uel’  Ue/te L?/tette/L? — JaV2ul” 

We may check that 8 is really non-dimensional. From the PDE, f must 
have the same dimensions as 0u/0t, i.e., [OT~']. The dimension of a is more 
intricate, but from the term auz, we know that uz, has dimensions [OL~?], 
and then a must have dimension [L?T~ '] to match the target [OT~']. In the 
expression for 6 we get [L?9T~'(L?T~'@)—1], which equals 1 as it should. 


B= 
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3.2.3 Jump boundary condition 


A classical one-dimensional heat conduction problem goes as follows. An 
insulated rod at some constant temperature Ug is suddenly heated from one 
end (x =0), modeled as a constant Dirichlet condition u(0,t) = U1 4 Uo at 
that end. That is, the boundary temperature jumps from Up to U; at t=0. All 
the other surfaces of the rod are insulated such that a one-dimensional model 
is appropriate, but we must explicitly demand uz(L,t) =0 to incorporate the 
insulation condition in the one-dimensional model at the end of the domain 
x = L. Heat cannot escape, and since we supply heat at x = 0, all of the 
material will eventually be warmed up to the temperature U1: u— Uj, as 
t> oo. 
The initial-boundary value problem reads 


Ou O7u 
u(x,0) = Uo, x € [0, L], (3.42) 
u(0,t) =U}, te (0,7), (3.43) 
2 ewe: tE(0,T]. (3.44) 
Ox 


The PDE (3.41) arises from the energy equation in solids and involves three 
physical parameters: the density o, the specific heat capacity parameter c, 
and the heat conduction coefficient (from Fourier’s law). Dividing by oc and 
introducing a = k/(oc) brings (3.41) on the standard form (3.36). We just 
use the a parameter in the following. 

The natural dimensionless temperature for this problem is 


since this choice makes ù € [0,1]. The reason is that u is bounded by the initial 
and boundary conditions (in the absence of a source term in the PDE), and 
we have u(Z,0) = 0, u(Z,0o) = 1, and u(0,t) = 1. 

The choice of te is as in the previous cases. We arrive at the dimensionless 
initial-boundary value problem 


aa ze (0,1), te (0,7), (3.45) 
u(z,0) =0, z€ [0,1], (3.46) 
u(0,t) =1, te (0,7, (3.47) 

2 u(1,t) =0, te (0,7). (3.48) 


da 
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The striking feature is that there are no physical parameters left in this 
problem. One simulation can be carried out for u(Z,t), see Figure 3.3, and the 
temperature in a rod of any material and any constant initial and boundary 
temperature can be retrieved by 


u(x,t) = Uo + (U1 — Up) a(a/L,ta/L?). 


0.0 0.2 0.4 0.6 0.8 1.0 


Fig. 3.3 Scaled temperature in an isolated rod suddenly heated from the end. 


3.2.4 Oscillating Dirichlet condition 


Now we address a heat equation problem where the temperature is oscillating 
on the boundary x = 0: 


Ou Oru 
BE aa? x €(0,L), te (0,7), (3.49) 
u(z,0) = Up, x € [0, L], (3.50) 
u(0,t) = Uo + Asin(wt), te (0,7), (3.51) 
LA =0, t € (0,T]. (3.52) 


Ox 
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One important physical application is temperature oscillations in the ground, 
either day and night variations at a short temporal and spatial scale, or 
seasonal variations in the Earth’s crust. An important modeling assumption 
is (3.52), which means that the boundary x = L is placed sufficiently far from 
x =0 such that the solution is much damped and basically constant so uz; = 0 
is a reasonable condition. 


Scaling issues. Since the boundary temperature is oscillating around the 
initial condition, we expect u € [Up — A, Uo + A]. The dimensionless temper- 
ature is therefore taken as 


such that u € [—1,1]. 

What is an appropriate time scale? There will be two time scales in- 
volved, the oscillations sin(wt) with period P = 27/w at the boundary and 
the “speed of diffusion”, or more specifically the “speed of heat conduction” 
in the present context, where te = x2/a is the appropriate scale, xe being 
the length scale. Choosing the right length scale is not obvious. As we shall 
see, the standard choice £e = L is not a good candidate, but to understand 
why, we need to examine the solution, either through simulations or through 
a closed-form formula. We are so lucky in this relatively simple pedagogical 
problem that one can find an exact solution of a related problem. 


Exact solution. As usual, investigating the exact solution of the model 
problem can illuminate the involved scales. For this particular initial-boundary 
value problem the exact solution as t + co (such that the initial condition 
u(x,0) = Uo is forgotten) and L > oo (such that (3.52) is certainly valid) can 
be shown to be 


w 
u(x,t) = Uo — Ae™™ sin(ba — wt), b= Va: (3.53) 
a 
This solution is of the form e~°* g(a — ct), i.e., a damped wave that moves 
to the right with velocity c and a damped amplitude e~>”. This is perhaps 
more easily seen if we make a rewrite 


u(x,t) = Up — Ae™®™ sin (b(a —ct)), c=w/b=V2aw, b=,/ A 


Time and length scales. The boundary oscillations lead to the time scale 
te = 1/w. The speed of the wave suggests another time scale: the time it 
takes to propagate through the domain, which is L/c, and hence te = L/c= 
L/V2aw. 

One can argue that L is not the appropriate length scale, because u is 
damped by e~°”. So, for x > 4/b, u is close to zero. We may instead use 
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1/b as length scale, which is the e-folding distance of the damping factor, 
and base te on the time it takes a signal to propagate one length scale, 
tz! = bc =w. Similarly, the time scale based on the “speed of diffusion” 
changes to tz! = b?a = jw if we employ 1/b as length scale. 

To summarize, we have three candidates for the time scale: te = L? /a 
(diffusion through the entire domain), te = 2/w (diffusion through a distance 
1/b where u is significantly different from zero), and te = 1/w (wave movement 
over a distance 1/0). 

Let us look at the dimensionless exact solution to see if it can help with 
the choice of scales. We introduce the dimensionless parameters 


P =bre = te x y= wte: 
V 2a 


The scaled solution becomes 


ulz, t; 8,7) =e °* sin(yt — 82). 
The three choices of y, implied by the three choices of te, are 


l, te=l1/w, 
Y=4 2, te=2/w, (3.54) 
282, te = L? /a, te= L 


The former two choices leave only 8 as parameter in U, and with ze = 1/b 
as length scale, 8 becomes unity, and there are no parameters in the dimen- 
sionless solution: 


u(Z,t) =e 7 sin(t— z7). (3.55) 


Therefore, £e = 1/b and te = 1/w (or te = 2/w, but the factor 2 is of no 
importance) are the most appropriate scales. 

To further argue why (3.55) demonstrates that these scales are preferred, 
think of w as large. Then the wave is damped over a short distance and there 
will be a thin boundary layer of temperature oscillations near x = 0 and 
little changes in u in the rest of the domain. The scaling (3.55) resolves this 
problem by using 1/b ~ w!/? as length scale, because then the boundary 
layer thickness is independent of w. The length of the domain can be chosen 
as, e.g., 4/b such that u +0 at the end « = L. The length scale 1/b helps us 
to zoom in on the part of u where significant changes take place. 

In the other limit, w small, b becomes small, and the wave is hardly damped 
in the domain (0, L] unless L is large enough. The imposed boundary condition 
on x = L in fact requires u to be approximately constant so its derivative 
vanishes, and this property can only be obtained if L is large enough to 
ensure that the wave becomes significantly damped. Therefore, the length 
scale is dictated by b, not L, and L should be adapted to b, typically L > 4/b if 
e~* ~ 0.018 is considered enough damping to consider ù ~ 0 for the boundary 
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condition. This means that x € [0,4/b] and then z € [0,4]. Increasing the 
spatial domain to [0,6] implies a damping e~° ~ 0.0025, if more accuracy is 
desired in the boundary condition. 


The scaled problem. Based on the discussion of scales above, we arrive at 
the following scaled initial-boundary value problem: 


u 2u 2 5 
ion, ze (0,4), F€(0,T], (3-56) 
u(z,0) = 0, z € [0,1], (3.57) 
u(0,t) = sin(t), te (0,7), (3.58) 
LA =0, te (0,T]. (3.59) 


Ox 
The coefficient in front of the second-derivative is 4 because 


tea ba 1 


1/2 w 2 


We may, of course, choose te = 2/w and get rid of the 5 factor, if desired, but 
then it turns up in (3.58) instead, as sin(2t). 

The boundary condition at 7 = L is only an approximation and relies on 
sufficient damping of ù to consider it constant (0/0Z = 0) in space. We could, 


therefore, assign the condition u=0 instead at z= L. 


Simulations. The file session. py contains a function solver_diffusion_FE 
for solving a diffusion equation in one dimension. This function can be used 
to solve the system (3.56)-(3.59), see diffusion_oscillatory_BC. 


Movie 3: Diffusion wave. https: //github.com/hplgit/scaling-book/raw/master/ 
doc/pub/book/htm1/mov-scaling/diffusion_osc_BC/movie.mp4 


3.3 Reaction-diffusion equations 


3.3.1 Fisher’s equation 


Fisher’s equation is essentially the logistic equation at each point for popula- 
tion dynamics (see Section 2.1.9) combined with spatial movement through 
ordinary diffusion: 


Ou 3u 
a Oe +ou(l1—u/M). (3.60) 


This PDE is also known as the KPP equation after Kolmogorov, Petrovsky, 
and Piskunov (who introduced the equation independently of Fisher). 
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Setting 
> x meee z u 
T= tSn US > 
Te te Ue 
results in 
Ou tea OP 
— = — + t.ou(1—-—u,.u/M). 


Balance of all terms. If all terms are equally important, the scales can be 
determined from demanding the coefficients to be unity. Reasoning as for the 
logistic ODE in Section 2.1.9, we may choose te = 1/0. Then the coefficient 
in the diffusion term dictates the length scale xe = V/t-a. A natural scale for 
u is M, since M is the upper limit of u in the model (cf. the logistic term). 
Summarizing, 


1 
ue = M, t=-, Te = sA 
Q Q 
and the scaled PDE becomes 
ðu uù 3 
DE E aR tu(l—u). (3.61) 


With this scaling, the length scale x, = \/a/g is not related to the domain 
size, so the scale is particularly relevant for infinite domains. 

An open question is whether the time scale should be based on the diffusion 
process rather than the initial exponential growth in the logistic term. The 
diffusion time scale means te = x2/a, but demanding the logistic term then 
to have a unit coefficient forces x20/a =1, which implies x, = \/a/o and 
te = 1/0. That is, equal balance of the three terms gives a unique choice of 
the time and length scale. 


Fixed length scale. Assume now that we fix the length scale to be L, either 
the domain size or some other naturally given length. With ze = L, te = 07}, 
Uc = M, we get 


ðù ouù 
-= u(l—u 3.62 
ap 7 Bags tUl), (3.62) 
where ( is a dimensionless number 
jE eat = 
oL? Ja 


The last equality demonstrates that 8 measures the ratio of the time scale for 
exponential growth in the beginning of the logistic process and the time scale 
of diffusion L?/a (i.e., the time it takes to transport a signal by diffusion 
through the domain). For small 6 we can neglect the diffusion and spatial 
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movements, and the PDE is essentially a logistic ODE at each point, while 
for large 8, diffusion dominates, and te should in that case be based on the 
diffusion time scale L?/a. This leads to the scaled PDE 


2 
ZA aga t6 u(L— 4), (3.63) 


showing that a large 6 encourages omission of the logistic term, because 
the point-wise growth takes place over long time intervals while diffusion is 
rapid. The effect of diffusion is then more prominent and it suffices to solve 
ug = Uzz. The observant reader will in this latter case notice that ue = M is an 
irrelevant scale for u, since logistic growth with its limit is not of importance, 
so we implicitly assume that another scale ue has been used, but that scale 
cancels anyway in the simplified PDE uj; = uaz. 


3.3.2 Nonlinear reaction-diffusion PDE 


A general, nonlinear reaction-diffusion equation in 1D looks like 


Ou 3u 

— =Q ; 3.64 

Y as yt flu) (3.64) 
By scaling the nonlinear reaction term f(u) as fef (uct), where fe is a charac- 
teristic size of f (u), typically the maximum value, one gets a non-dimensional 
PDE like 


Ou tean tefes, _ 
Ot a2 OB? uc F{uett). 


The characteristic size of u can often be derived from boundary or initial 
conditions, so we first assume that ue is given. This fact uniquely determines 
the space and time scales by demanding that all three terms are equally 
important and of unit size: 


E Ra 
aa as fe 

The corresponding PDE reads 
ða Pu - 
an fluu). (3.65) 


If xe is based on some known length scale L, balance of all three terms 
can be used to determine ue and te: 


L Lf 
t= Ue = : 
a a 
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This scaling only works if f is nonlinear, otherwise ue cancels and there is no 
freedom to constrain this scale. 

With given L and uc, there are two choices of te since it can be based on 
the diffusion or the reaction time scales. With the reaction scale, te = uc/ fe, 
one arrives a the PDE 

du u - 


BE = Baa t f (ucu), (3.66) 


where 


B= Auc _ Ue /fe 
Lfe L?/a 
is a dimensionless number reflecting the ratio of the reaction time scale and 


the diffusion time scale. On the contrary, with the diffusion time scale, te = 
L? /a, the scaled PDE becomes 


2 oy 
ap ga tt BO Flucti). (3.67) 


The size of 8 in an application will determine which of the scalings that is 
most appropriate. 


3.4 The convection-diffusion equation 


3.4.1 Convection-diffusion without a force term 


We now add a convection term v- Vu to the diffusion equation to obtain the 
well-known convection-diffusion equation: 


—+vy-Vu=aV7u, 2,y,2€2, t€ (0,7). (3.68) 


The velocity field v is prescribed, and its characteristic size V is normally 
clear from the problem description. In the sketch below, we have some given 
flow over a bump, and u may be the concentration of some substance in the 
fluid. Here, V is typically max, u(y). The characteristic length L could be 
the entire domain, L = c+, or the height of the bump, L = D. (The latter 
is the important length scale for the flow.) 
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Inserting 
a Do a Y . Z F t > v _ u 
x F Z ? U > U 
Te Y Ye Ze te V Uc 
in (3.68) yields 
Uc OU UcV Quc = 


ö-Vü= 


Eon : —V7u, 2,9,2€ Q, te (0,T). 
For ue we simply introduce the symbol U, which we may estimate from an 
initial condition. It is not critical here, since it vanishes from the scaled equa- 
tion anyway, as long as there is no source term present. With some velocity 
measure V and length measure L, it is tempting to just let te = L/V. This is 
the characteristic time it takes to transport a signal by convection through 
the domain. The alternative is to use the diffusion length scale te = L? /a. A 
common physical scenario in convection-diffusion problems is that the con- 
vection term v- Vu dominates over the diffusion term aV?u. Therefore, the 
time scale for convection (L/V) is most appropriate of the two. Only when the 
diffusion term is very much larger than the convection term (corresponding 
to very small Peclet numbers, see below) te = L? /a is the right time scale. 
The non-dimensional form of the PDE with te = L/V becomes 


Ot 4 0-Va= Pe Va, z,4,2 € 2, te (0,7), (3.69) 
where Pe is the Peclet number, 
L 
Pe= ir 3 
a 


Estimating the size of the convection term v: Vu as VU/L and the diffusion 
term aV7u as aU/L?, we see that the Peclet number measures the ratio of 
the convection and the diffusion terms: 


_ convection VU/L LV 
diffusion = aU/L?2 a” 
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In case we use the diffusion time scale te = L?/a, we get the non- 
dimensional PDE 


—+Pev-Vi=V*t, £,9,2€ 2, t€ (0,T]. (3.70) 


Discussion of scales and balance of terms in the PDE 


We see that (3.69) and (3.70) are not equal, and they are based on two 
different time scales. For moderate Peclet numbers around 1, all terms 
have the same size in (3.69), i.e., a size around unity. For large Peclet 
numbers, (3.69) expresses a balance between the time derivative term 
and the convection term, both of size unity, and then there is a very 
small Pe~'V?u term because Pe is large and Vu should be of size 
unity. That the convection term dominates over the diffusion term is 
consistent with the time scale te = L/V based on convection transport. 
In this case, we can neglect the diffusion term as Pe goes to infinity and 
work with a pure convection (or advection) equation 
Ou 


ate Va=0. 


For small Peclet numbers, Pe~'V?% becomes very large and can only 


be balanced by two terms that are supposed to be unity of size. The 
time-derivative and/or the convection term must be much larger than 
unity, but that means we use suboptimal scales, since right scales imply 
that 0u/Ot and v: Vü are of order unity. Switching to a time scale based 
on diffusion as the dominating physical effect gives (3.70). For very small 
Peclet numbers this equation tells that the time-derivative balances 
the diffusion. The convection term 0- Vu is around unity in size, but 
multiplied by a very small coefficient Pe, so this term is negligible in 
the PDE. An approximate PDE for small Peclet numbers is therefore 


Scaling can, with the above type of reasoning, be used to neglect 
terms from a differential equation under precise mathematical condi- 
tions. 
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3.4.2 Stationary PDE 


Suppose the problem is stationary and that there is no need for any time 
scale. How is this type of convection-diffusion problem scaled? We get 


VU _ =_ aU — ee 
ore a 
or 
o-Vu=Pe!V7u. (3.71) 


This scaling only “works” for moderate Peclet numbers. For very small or 
very large Pe, either the convection term 0- Vu or the diffusion term V?u 
must deviate significantly from unity. 

Consider the following 1D example to illustrate the point: v = vi, v > 0 
constant, a domain [0, L], with boundary conditions u(0) =0 and u(L) = UŁ. 
(The vector į is a unit vector in x direction.) The problem with dimensions 
is now 


d d 
ieee gga? ZEI), a(0)=0, a(1)=1, 
if we choose U = Uz. The solution of the scaled problem is 
ae 1 _ etPe 
aa ae 


Figure 3.4 indicates how u depends on Pe: small Pe values give approximately 
a straight line while large Pe values lead to a boundary layer close to x = 1, 
where the solution changes very rapidly. 

We realize that for large Pe, 


which are consistent results with the PDE, since the double derivative term 
is multiplied by Pe‘. For small Pe, 
max du xl, max d'n x 
a dt ° z dz? 
which is also consistent with the PDE, since an almost vanishing second-order 
derivative is multiplied by a very large coefficient Pe~!. However, we have a 
problem with very large derivatives of u when Pe is large. 


0, 
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Fig. 3.4 Solution of scaled problem for 1D convection-diffusion. 


To arrive at a proper scaling for large Peclet numbers, we need to remove 
the Pe coefficient from the differential equation. There are only two scales at 
our disposals: ue and ze for u and 2, respectively. The natural value for ue is 
the boundary value Uz; at «= L. The scaling of Vugz = quzz then results in 


duù œa u = 
di Va, dz?’ 


where L= L/2¢. Choosing the coefficient a/(Va-) to be unity results in the 
scale £e =a/V, and L becomes Pe. The final, scaled boundary-value problem 
is now 


di Pa 
== 55 FE(0,Pe), w0)=0, (Pe) =1, 
with solution 
aa 1—e* 
ult) = Pe" 


Figure 3.5 displays u for some Peclet numbers, and we see that the shape 
of the graphs are the same with this scaling. For large Peclet numbers we 
realize that ū and its derivatives are around unity (1—eP® ~ —eP®), but for 
small Peclet numbers dū/dz ~ Pe~*. 
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Fig. 3.5 Solution of scaled problem where the length scale depends on the Peclet 
number. 


The conclusion is that for small Peclet numbers, £e = L is an appropriate 
length scale. The scaled equation Pew’ = u” indicates that u” ~ 0, and the 
solution is close to a straight line. For large Pe values, xe = a/V is an appro- 
priate length scale, and the scaled equation U’ = U” expresses that the terms 
w and u” are equal and of size around unity. 


3.4.3 Convection-diffusion with a source term 


Let us add a force term f(x,t) to the convection-diffusion equation: 


o 
S to Vu= aV2ut f. (3.72) 
The scaled version reads 
dù tV_ = teasa tefez 
ae .Va= Vv fede 
a e 


We can base te on convective transport: te = L/V. Now, uc could be chosen 
to make the coefficient in the source term unity: ue = tefe = Lfe/V. This 
leaves us with 
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aðu _ _ = 
OY 4 b-Va= Pe !V2a4 fF. 
ot 
In the diffusion limit, we base te on the diffusion time scale: te = L?/a, 
and the coefficient of the source term set to unity determines ue according to 


L*f. L 
fe si = wui = fe : 
Ue a 


The corresponding PDE reads 


ot pa Vut f, 
Ot 
so for small Peclet numbers, which we have, the convective term can be 
neglected and we get a pure diffusion equation with a source term. 
What if the problem is stationary? Then there is no time scale and we get 


Vue = Uc ao_ = 
1 oye T Vut fef, 
or 
é, 2 -z 
ovas Pe va + E, 
Vue 


Again, choosing ue such that the source term coefficient is unity leads to 
Uc = feL/V. Alternatively, uc can be based on the initial condition, with 
similar results as found in the sections on the wave and diffusion PDEs. 
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Chapter 4 


Advanced partial differential 
equation models 


This final chapter addresses more complicated PDE models, including linear 
elasticity, viscous flow, heat transfer, porous media flow, gas dynamics, and 
electrophysiology. A range of classical dimensionless numbers are discussed 
in terms of the scaling. 


4.1 The equations of linear elasticity 


To the best of the authors’ knowledge, it seems that mathematical models in 
elasticity and structural analysis are almost never non-dimensionalized. This 
is probably due to tradition, but the following sections will demonstrate the 
usefulness of scaling also in this scientific field. 

We start out with the general, time-dependent elasticity PDE with variable 
material properties. Analysis based on scaling is used to determine under 
what circumstances the acceleration term can be neglected and we end up 
with the widely used stationary elasticity PDE. Scaling of different types 
of boundary conditions is also treated. At the end, we scale the equations 
of coupled thermo-elasticity. All the models make the assumption of small 
displacement gradients and Hooke’s generalized constitutive law such that 
linear elasticity theory applies. 


4.1.1 The general time-dependent elasticity problem 


The following vector PDE governs deformation and stress in purely elastic 
materials, under the assumption of small displacement gradients: 


2u 
0B = V((A+u)Y u) +V- (uVu)+ of. (4.1) 
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Here, u is the displacement vector, o is the density of the material, À and u 
are the Lame elasticity parameters, and f is a body force (gravity, centrifugal 
force, or similar). 

We introduce dimensionless variables: 


= =i De aR ~ Y . 4 ~_f 
= t ; 
U=U. U, T T Yy L A T’ i 


Also the elasticity parameters and the density can be scaled, if they are not 
constants, 


z _À `_H- __ 2 
À a, aoe L=—_, = —, 
c He Qc 
where the characteristic quantities are typically spatial maximum values of 


the functions: 


Ace =mMaxA, Pe=Maxp, Oc =mMaxga. 
T,Y,Z T,Y,Z T,Y,Z 


Finally, we scale f too (if not constant): 


f=fo'f, fe= max ||fI]. 
aA EA 


Inserting the dimensionless quantities in the governing vector PDE results 
in 


Pu te 7 a tHe i te Ee fie 
IP 7 22 V(AV +d) + ag V lEn (Vù) + aD f. (4.2) 


We may choose te to make the coefficient in front of any of the spatial deriva- 
tive terms equal unity. Here we choose the u term, which implies 


aee, 
Hec 


The scale for u can be chosen from an initial displacement or by making the 
coefficient in front of the f term unity. The latter means 


Uc = uz efb? . 


As discussed later, in Section 4.1.4, this might not be the desired uc in 
applications. 
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The resulting dimensionless PDE becomes 


DOT Ape nn cate, oie 
Oy = V((GAFM)V- 0) +V- (Vu) + af. (4.3) 
The only dimensionless parameter is 
yee 
He 


If the source term is absent, we must use the initial condition or a known 
boundary displacement to determine ue. 


Software. Given software for (4.1), we can simulate the dimensionless prob- 
lem by setting o = 0, A= GBA, and u= A. 


4.1.2 Dimensionless stress tensor 


The stress tensor ø is a key quantity in elasticity and is given by 


o =AV-ul+p(Vut(Vu)?). 


This ø can be computed as soon as the PDE problem for u has been solved. 
Inserting dimensionless variables on the right-hand side of the above relation 
gives 


o = etic LV -U+ pec Lt fi(Va + (Vu)*) 
Sinh. (sav a+ (Vat (Vu)")) 
The coefficient on the right-hand side, peucLT!, has dimension of stress, 
since (according to the second table in Section 1.1.2) [MT~?L~1)(L)(L~)] = 


[MT~?L~ 14], which is the dimension of stress. The quantity ycuc L+ is there- 
fore the natural scale of the stress tensor: 


-1 
» Co=MctcL, 
and we have the dimensionless stress-displacement relation 


o=BA\V-ut+p(Vut(Vu)"). (4.4) 
4.1.3 When can the acceleration term be neglected? 


A lot of applications of the elasticity equation involve static or quasi-static 
deformations where the acceleration term ou;; is neglected. Now we shall see 
under which conditions the quasi-static approximation holds. 
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The further discussion will need to look into the time scales of elastic 
waves, because it turns out that the chosen te above is closely linked to the 
propagation speed of elastic waves in a homogeneous body without body 
forces. A relevant model for such waves has constant o, A, and u, and no 
force term: 


2 
oes =(A+p)VV-utpV7u. (4.5) 
S waves. Let us take the curl of this PDE and notice that the curl of a 
gradient vanishes. The result is 


2 
ap Y xu= VV x U, 


i.e., a wave equation for V x u. The wave velocity is 


_ JE 
cg=,/-. 
Q 


The corresponding waves are called S waves!. The curl of a displacement field 
is closely related to rotation of continuum elements. S waves are therefore 
rotation waves, also sometimes referred to as shear waves. 

The divergence of a displacement field can be interpreted as the volume 
change of continuum elements. Suppose this volume change vanishes, V -u = 
0, which means that the material is incompressible. The elasticity equation 
then simplifies to 


Pu 

ap =° 
so each component of the displacement field in this case also propagates as 
a wave with speed ch. The time it takes for such a wave to travel one charac- 
teristic length L is L/cg, i.e., Ly/o/p , which is nothing but our characteristic 
time te. 


207 u, 


P waves. We may take the divergence of the PDE instead and notice that 
V-V = V? so 


82 
a un cbV?V.-u, 


A+ 2p 
cp= 
Q 


lhttps://en.wikipedia.org/wiki/S-wave 


with wave velocity 
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That is, the volume change (expansion/compression) propagates as a wave 
with speed cp. These types of waves are called P waves”. Other names are 
pressure and expansion/compression waves. 

Suppose now that V x u=0, i.e., there is no rotation (“shear”) of contin- 
uum elements. Mathematically this condition implies that V?u = V(V-u) 
(see any book on vector calculus or Wikipedia®). Our model equation (4.5) 
then reduces to 


Oru 

Ot 
which is nothing but a wave equation for the expansion component of the 
displacement field, just as (4.1.3) is for the shear component. 


= c3,V*u, 


Time-varying load. Suppose we have some time-varying boundary con- 
dition on u or the stress vector (traction), with a time scale 1/w (some 
oscillating movement that goes like sinwt, for instance). We choose te = 1/w. 
The scaling now leads to 


OO Wy, sea tee Sa le ae! 20 
Yge = VUBAT MV u)+V- (VU) + of. 


where we have set 


Uc = uzt feL’ oc, 


as before, and y is a new dimensionless number, 


i oL’ w? E Ly/ 0c/ he ? 
q Le 1/w 
The last rewrite shows that ,/y is the ratio of the time scale for S waves and 
the time scale for the forced movement on the boundary. The acceleration 
term can therefore be neglected when y < 1, i.e., when the time scale for 
movement on the boundary is much larger than the time it takes for the S 
waves to travel through the domain. Since the velocity of S waves in solids is 


very large and the time scale correspondingly small, y <1 is very often the 
case in applications involving structural analysis. 


4.1.4 The stationary elasticity problem 


Scaling of the PDE. We now look at the stationary version of (4.1) where 
the ou: term is removed. The first step in the scaling is just inserting the 
dimensionless variables: 


*https://en.wikipedia. org/wiki/P-wave 
Shttps://en.wikipedia. org/wiki/Vector_calculus_identities 
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0 = Lo ucV((AcA + Mel) Vt) + L7 ’?ucheV (VU) + Oc feof - 
Dividing by Luche gives 
bose ee SN a = Totes 
0=V((86A+0)V -u)+V (Vu) + ——of. 
Uche 
Choosing uc = oL? fe/He leads to 
V((BA+f)V-a)+V-(iVu)+ of =0. (4.6) 
A homogeneous material with constant A, 4, and @ is an interesting case 
(this corresponds to He = H, Ac = À, Oc = 0, O=A=pP=1): 


(1+ 6)V(V-a)+V7ua)+f =0. (4.7) 


pees ee: 
H Cs 


It shows that in standard, stationary elasticity, å/u is the only significant 
physical parameter. 


Now £ is defined as 


Remark on the characteristic displacement. 


The presented scaling may not be valid for problems where the geometry 
involves some dimensions that are much smaller than others, such as for 
beams, shells, or plates. Then one more length scale must be defined which 
gives us non-dimensional geometrical numbers. Global balances of moments 
and loads then determine how characteristic displacements depend on these 
numbers. As an example, consider a cantilever beam of length L and square- 
shaped cross section of width W, deformed under its own weight. From beam 
theory one can derive uc = $0gL75*/E, where 6 = L/W (g is the acceleration 
of gravity). If we consider E to be of the same size as A, this implies that 
y~ 672. So, it may be wise to prescribe a ue in elasticity problems, perhaps 
from formulas as shown, and keep y in the PDE. 


Scaling of displacement boundary conditions. A typical boundary con- 
dition on some part of the boundary is a prescribed displacement. For sim- 
plicity, we set u = Up for a constant vector Ug as boundary condition. With 
Uc = oL? fe/ u, we get the dimensionless condition 

Uy _ HU 

Uc oL? fe 

In the absence of body forces, the expression for uc has no meaning (fe = 0), 


so then uc = |Uo] is a better choice. This gives the dimensionless boundary 
condition 


Sı 


Uo 
|Uo|’ 


which is the unit vector in the direction of Up. The new ue changes the 


u = 
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coefficient in front of the body force term, if that term is present, to the 
dimensionless number 
_ L of. G 

L|Up| 


Scaling of traction boundary conditions. The other type of common 
boundary condition in elasticity is a prescribed traction (stress vector) on 
some part of the boundary: 


ô 


o:n= T, 


where, to make it simple, we take Jo as a constant vector. From Section 4.1.2 
we have a stress scale oc = puc/L, but we may alternatively use |To| as stress 
scale. In that case, 
_ Fo 

[Tol 


o'n 


which is a unit vector in the direction of Tg. Many applications involve large 
traction free areas on the boundary, on which we simply have o-n =0. 


4.1.5 Quasi-static thermo-elasticity 


Heating solids gives rise to expansion, i.e., strains, which may cause stress 
if displacements are constrained. The time scale of temperature changes are 
usually much larger than the time scales of elastic waves, so the stationary 
equations of elasticity can be used, but a term depends on the temperature, 
so the equations must be coupled to a PDE for heat transfer in solids. The 
resulting system of PDEs is known as the equations of thermo-elasticity and 
reads 


V((A+4)V -u)+V-(uVu)=aVT-—-of, (4.8) 
ge =V- (KYT) + ofr, (4.9) 


where T is the temperature, a is a coefficient of thermal expansion, c is a heat 
capacity, k is the heat conduction coefficient, and fr is some heat source. 
The density o is strictly speaking a function of T and the stress state, but 
a widely used approximation is to consider ọ as a constant. Most thermo- 
elasticity applications have fr = 0, so we drop this term. Most applications 
also involve some heating from a temperature level Tọ to some level Ty + AT. 
A suitable scaling for T is therefore 


T-T 
AT ’ 
so that T € [0,1]. The elasticity equation has already been scaled and so has 


the diffusion equation for T. We base the time scale on the diffusion, i.e., the 
thermal conduction process: 


T= 


t= ocL? fre. 
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We imagine that « is scaled as K = k/K . The dimensionless PDE system then 
becomes 

=VT—eðf, (4.10) 
- =V- (KVT). (4.11) 
Here we have chosen ue such that the “heating source term” has a unit 


coefficient, acknowledging that this thermal expansion balances the stress 
terms with u. The corresponding displacement scale is 


aLAT 
Ue = : 
Hec 
The dimensionless number in the body force term is therefore 
pa Leche 
aAT’ 


which measures the ratio of the body force term and the “heating source 
term”. 

A homogeneous body with constant o, A, u, c, and « is common. The PDE 
system reduces in this case to 


V((1+ B)V-a)+V2a) = VT —-ef, (4.12) 
OT Spe 
apa We: (4.13) 


In the absence of body forces, 8 is again the key parameter. 

The boundary conditions for thermo-elasticity consist of the conditions 
for elasticity and the conditions for diffusion. Scaling of such conditions are 
discussed in Section 3.2 and 4.1.4. 


4.2 The Navier-Stokes equations 


This section shows how to scale various versions of the equations governing 
incompressible viscous fluid flow. We start with the plain Navier-Stokes equa- 
tions without body forces and progress with adding the gravity force and a 
free surface. We also look at scaling low Reynolds number flow and oscillating 
flows. 
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4.2.1 The momentum equation without body forces 


The Navier-Stokes equations for incompressible viscous fluid flow, without 
body forces, take the form 


o( Ge twvu) =—Vp+pV2u, (4.14) 


V-u=0. (4.15) 


The primary unknowns are the velocity u and the pressure p. Moreover, o is 
the fluid density, and ys is the dynamic viscosity. 


Scaling. We start, as usual, by introducing a notation for dimensionless 
independent and dependent variables: 


£ y z t u p 


L >) y L ’ L 9 tz 3 uc 9 De Fi 
where L is some characteristic distance, te is some characteristic time, ue is 
a characteristic velocity, while pe is a characteristic pressure. Inserted in the 
equations, 


Uc OU u2 _ = Pea Uc =2- 
ahs Cie =- pain 4.1 
($ a t Tu vu) 7 VP+ paeV ù, (4.16) 
ZV u=0. (4.17) 


For the velocity it is common to just introduce some U for ue. This U is 
normally implied by the problem description. For example, in the flow con- 
figuration below, with flow over a bump, we have some incoming flow with 
a profile u(y) and U can typically be chosen as U = max, v(y). The height 
of the bump influences the wake behind the bump, and is the length scale 
that really impacts the flow, so it is natural to set L = D. For numerical 
simulations in a domain of finite extent, [0,c+ 4, c must be large enough 
to avoid feedback on the inlet profile, and £ must be large enough for the 
type of outflow boundary condition used. Ideally, c,?— oo, so none of these 
parameters are useful as length scales. 
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For flow in a channel or tube, we also have some inlet profile, e.g., v(r) 
in a tube, where r is the radial coordinate. A natural choice of characteristic 
velocity is U = v(0) or to let U be the average flow, i.e., 


1 R 
U= =a | 2rv(r)rdr, 


if R is the radius of the tube. Other examples may be flow around a body, 
where there is some distant constant inlet flow u =  Upo?, for instance, and 
U = Up is an obvious choice. We therefore assume that the flow problem 
itself brings a natural candidate for U. 

Having a characteristic distance L and velocity U, an obvious time measure 
is L/U so we set te = L/U. Dividing by the coefficient in front of the time 
derivative term, creates a pressure term 


The coefficient suggest a choice pe = oU? if the pressure gradient term is 
to have the same size as the acceleration terms. This pe is a very common 
pressure scale in fluid mechanics, arising from Bernoulli’s equation 
1 
p+ zeuu = const 
for stationary flow. 


Dimensionless PDEs and the Reynolds number. The discussions so far 
result in the following dimensionless form of (4.14) and (4.15): 


—+u-Vu=—Vp+Re tV’ u, (4.18) 
V-a=0, (4.19) 
where Re is the famous Reynolds number, 


_ wb _UL 
oa 


Re 
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The latter expression makes use of the kinematic viscosity v = y/o. For vis- 
cous fluid flows without body forces there is hence only one dimensionless 
number, Re. 

The Reynolds number can be interpreted as the ratio of convection and 
viscosity: 


convection |ou-Vu| @U?/L UL R 
= = =Re. 


viscosity  |uV?2ul a pU/L2 v 
(We have here used that Vu goes like U/L and V?u goes like U/L?.) 


4.2.2 Scaling of time for low Reynolds numbers 


As we discussed in Section 3.4 for the convection-diffusion equation, there is 
not just one scaling that fits all problems. Above, we used te = L/U, which 
is appropriate if convection is a dominating physical effect. In case the con- 
vection term ou: Vu is much smaller than the viscosity term uVĉ?u, i.e., 
the Reynolds number is small, the viscosity term is dominating. However, if 
the scaling is right, the other terms are of order unity, and Re’ 'V?u must 
then also be of unit size. This fact implies that V?ū must be small, but then 
the scaling is not right (since a right scaling will lead to ù and its derivatives 
around unity). Such reasoning around inconsistent size of terms clearly points 
to the need for other scales. 

In the low-Reynolds number regime, the diffusion effect of V7u is dominat- 
ing, and we should use a time scale based on diffusion rather than convection. 
Such a time scale is te = L? /(u/o) = L? /v. With this time scale, the dimen- 
sionless Navier-Stokes equations look like 


au E ae 
op + Rea: Va =—Vp+ Va, (4.20) 


V-u=0. (4.21) 


As stated in the box in Section 3.4, (4.20) is the appropriate PDE for very 
low Reynolds number flow and suggests neglecting the convection term. If 
the flow is also steady, the time derivative term can be neglected, and we end 
up with the so-called Stokes problem for steady, slow, viscous flow: 


This flow regime is also known as Stokes’ flow or creeping flow. 
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4.2.3 Shear stress as pressure scale 


Instead of using the kinetic energy oU? as pressure scale, one can use the 
shear stress U/L (U/L reflects the spatial derivative of the velocity, which 
enters the shear stress expression wOu/Oy). Using U as velocity scale, L/U 
as time scale, and U/L as pressure scale, results in 


Re (2z +a- va) =-Vp+V’u. (4.24) 


Low Reynolds number flow now suggests neglecting both acceleration terms. 


4.2.4 Gravity force and the Froude number 
We now add a gravity force to the momentum equation (4.14): 
ðu 2 
o gte V" =—-Vp+uV~u-— ogk, (4.25) 


where g is the acceleration of gravity, and k is a unit vector in the oppo- 
site direction of gravity. The new term takes the following form after non- 
dimensionalization: 


mot > UP 
where Fr is the dimensionless Froude number, 
U 
vLg 


This quantity reflects the ratio of inertia and gravity forces: 


Fr= 


ju-Vul oU?/L 


Fr?. 
|og| og 


4.2.5 Oscillating boundary conditions and the Strouhal 
number 


Many flows have an oscillating nature, often arising from some oscillating 
boundary condition. Suppose such a condition, at some boundary x = const, 
takes the specific form 


u = Usin(wt)t. 
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The dimensionless counterpart becomes 


L 
Uu =U sin(w—t)e, 
Wgt 
if te = L/U is the appropriate time scale. This condition can be written 


u = sin(Stt), (4.26) 
where St is the Strouhal number, 


wL 
Se 
The two important dimensionless parameters in oscillating flows are then the 
Reynolds and Strouhal numbers. 

Even if the boundary conditions are of steady type, as for flow around a 
sphere or cylinder, the flow may at certain Reynolds numbers get unsteady 
and oscillating. For 10? < Re < 10’, steady inflow towards a cylinder will 
cause vortex shedding: an array of vortices are periodically shedded from the 
cylinder, producing an oscillating flow pattern and force on the cylinder. The 
Strouhal number is used to characterize the frequency of oscillations. The 
phenomenon, known as von Karman vortex street, is particularly important 
if the frequency of the force on the cylinder hits the free vibration frequency of 
the cylinder such that resonance occurs. The result can be large displacements 
of the cylinder and structural failure. A famous case in engineering is the 
failure of the Tacoma Narrows suspension bridge* in 1940, when wind-induced 
vortex shedding caused resonance with the free torsional vibrations of the 
bridge. 


St (4.27) 


4.2.6 Cavitation and the Euler number 


The dimensionless pressure in (4.18) made use of the pressure scale pe = oU ?. 
This is an appropriate scale if the pressure level is not of importance, which is 
very often the case since only the pressure gradient enters the flow equation 
and drives the flow. However, there are circumstances where the pressure 
level is of importance. For example, in some flows the pressure may become 
so low that the vapor pressure of the liquid is reached and that vapor cavities 
form (a phenomenon known as cavitation). A more appropriate pressure scale 
is then pe = Poo — Py, Where Poo is a characteristic pressure level far from 
vapor cavities and p, is the vapor pressure. The coefficient in front of the 
dimensionless pressure gradient is then 


fhttps ://en.wikipedia.org/wiki/Tacoma_Narrows_Bridge_(1940) 
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Poo — Pu 
oU? ` 
Inspired by Bernoulli’s equation p+ sou -u = const in fluid mechanics, a fac- 
tor 5 is often inserted in the denominator. The corresponding dimensionless 
number, 


_ Poo — Pv 


Eu= 2 tv 
30U? 


(4.28) 


is called the Euler number. The pressure gradient term now reads ZEuVp. 
The Euler number expresses the ratio of pressure differences and the kinetic 
energy of the flow. 


4.2.7 Free surface conditions and the Weber number 


At a free surface, z = 7(a,y,t), the boundary conditions are 


On 
= . 4.29 
w =n +u: Vn, ( ) 
an an 


where w is the velocity component in the z direction, po is the atmospheric 
air pressure at the surface, and o represents the surface tension. The approx- 
imation in (4.30) is valid under small deformations of the surface. 

The dimensionless form of these conditions starts with inserting the di- 
mensionless quantities in the equations: 


=) Lq ag 
ud = aE tUe: Vi, 
aal an | On 
PeP X-T? \ ax?" ay) 


The characteristic length L is usually taken as the depth of the fluid when 
the surface is flat. We have used p = (p -— po)/pe for making the pressure 
dimensionless. Using ue = U, te = L/U, and pe = oU?, results in 


w= ay te Vi (4.31) 


25 825 
px —We* (Satie) (4.32) 
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where We is the Weber number, 


U2L 
W=, (4.33) 


oO 


The Weber number measures the importance of surface tension effects and 
is the ratio of the pressure scale oU? and the surface tension force per area, 
typically o/R, in a 2D problem, which has size g/L. 


4.3 Thermal convection 


Temperature differences in fluid flow cause density differences, and since cold 
fluid is heavier than hot fluid, the gravity force will induce flow due to den- 
sity differences. This effect is called free thermal convection and is key to our 
weather and heating of rooms. Forced convection refers to the case where 
there is no feedback from the temperature field to the motion, i.e., tempera- 
ture differences do not create motion. This fact decouples the energy equation 
from the mass and momentum equations. 


4.3.1 Forced convection 


The model governing forced convection consists of the Navier-Stokes equa- 
tions and the energy equation for the temperature: 


a 
o (Z +u- vu) = —Vp+pV2u-—ogk, (4.34) 
V-u=0, (4.35) 
T 
0 (5 +u-V7) =KV7T. (4.36) 


The symbol T is the temperature, c is a heat capacity, and « is the heat 
conduction coefficient for the fluid. The PDE system applies primarily for 
liquids. For gases one may need a term —pV -u for the pressure work in 
(4.36) as well as a modified equation of continuity (4.35). 

Despite the fact that o depends on T, we treat o as a constant 09. The ma- 
jor effect of the o(T) dependence is through the buoyancy effect caused by the 
gravity term —0(T)gk. It is common to drop this term in forced convection, 
and assume the momentum and continuity equations to be independent of 
the temperature. The flow is driven by boundary conditions (rather than den- 
sity variations as in free convection), from which we can find a characteristic 
velocity U. 


114 4 Advanced partial differential equation models 


Dimensionless parameters are introduced as follows: 


Other coordinates are also scaled by L. The characteristic temperature Te is 
chosen as some range AT’, which depends on the problem and is often given by 
the thermal initial and/or boundary conditions. The reference temperature 
To is also implied by prescribed conditions. Inserted in the equations, we get 


: TU OT UT. a. oq KT. = 
a ar ae L 


Making each term in each equation dimensionless reduces the system to 


Ou - - = 

ap ti Va= _Vp+Re!V7a, (4.37) 
V-u=0, (4.38) 

Sta OF = Pe WF. (4.39) 


The two dimensionless numbers in this system are given by 


UL UL 
Pe= VZ Re= ge 


KE v oo 

The Peclet number is here defined as the ratio of the convection term for 
heat oọọcU AT/L and the heat conduction term KU/L?. The fraction «/(oo¢) 
is known as the thermal diffusivity, and if this quantity is given a symbol a, 
we realize the relation to the Peclet number defined in Section 3.4. 


4.3.2 Free convection 


Governing equations. The mathematical model for free thermal convec- 
tion consists of the Navier-Stokes equations coupled to an energy equation 
governing the temperature: 
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ð 
0 (Sr +u- Vu) = —Vp+ pV?u— ogk, (4.40) 
Op 
pr +V (gu) =0, (4.41) 
OT 2 
go( = +u- VT ) =KV?T + 2ueizeis, (4.42) 


where Einstein’s summation convention is implied for the e,je;; term. The 
symbol T is the temperature, c is a heat capacity, « is the heat conduction 
coefficient for the fluid. In free convection, the gravity term —0(T)gk is es- 
sential since the flow is driven by temperature differences and the fact that 
hot fluid rises while cold fluid falls. 

For slightly compressible gas flow a term —pV -u may be needed in (4.42). 


Heating by viscous effects. We have also included heating of the fluid 
due to the work of viscous forces, represented by the term 2uEijEij, where 
Eij is the strain-rate tensor in the flow, defined by 


hes 1 Ou; r Ou; = 1 ; T 
ot = 2 (= i iy 4 g (Vu (Vu) ) 

where u; is the velocity in direction of x; (i = 1,2,3 measures the space 
directions). The term 2ye;;¢;; is actually much more relevant for forced con- 
vection, but was left out in Section 4.3.2 for mathematical simplicity. Heating 
by the work of viscous forces is often a very small effect and can be neglected, 
although it plays a major role in forging and extrusion of metals where the 
viscosity is very large (such processes require large external forces to drive 
the flow). The reason for including the work by viscous forces under the 
heading of free convection, is that we want to scale a more complete, gen- 
eral mathematical model for mixed force and free convection, and arrive at 
dimensionless numbers that can tell if this extra term is important or not. 


Relation between density and temperature. Equations (4.40) has already 
been made dimensionless in the previous section. The major difference is now 
that o is no longer a constant, but a function of T. The relationship between 
o and T is often taken as linear, 


0 = 00 — 20b (T — To), 


where 


is known as the thermal expansion coefficient of the fluid, and ọọ is a reference 
density when the temperature is at Tọ. 
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The Boussinesq approximation. A very common approximation, called 
the Boussinesq approximation, is to neglect the density variations in all terms 
except the gravity term. This is a good approximation unless the change in 
o is large. With the linear o(T) formula and the Boussinesq approximation, 
(4.40)-(4.42) take the form 


Ou 
00 (= +u- vu) = —Vp+pV*u— (00 — 008(T — To))gk, (4.43) 
OT 2 
0c oe +u: VT | = KVT +H Qn ei; - (4.45) 


A good justification of the Boussinesq approximation is provided by Tritton 
(9, Ch. 13). 


Scaling. Dimensionless variables are introduced as 


-n g L ` u p = T-To 
Tora ue a 
The dimensionless y and z coordinates also make use of L as scale. As in 
forced convection, we assume the characteristic temperature level To and the 
scale AT are given by thermal boundary and/or initial conditions. Contrary 
to Sections 4.2 and 4.3.2, U is now not given by the problem description, but 
implied by AT. 
Replacing quantities with dimensions by their dimensionless counterparts 
results in 


U? ðu V? -_ PU ae £ 
— Pevp+ Ta V7u— oogk + oo8T eT gk, 


TU OT UTe_ == klesa t aU 


These equations reduce to 


ap ta Va= —Vp+Re 'V?u—Fr7*k+4Tk, (4.46) 
V-u=0, (4.47) 

OT = aig: pea —]1 = 2 Tr sane =] 

ar a a V T + 2ôēijEij « (4.48) 


The dimensionless numbers, in addition to Re and Fr, are 
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gBLAT Pet _* j= HU 
U2 ’ ~ ogcU L’ ~ LogcAT ` 


The y number measures the ratio of thermal buoyancy and the convection 
term: 
_ 9098AT _ gBLAT 
ooU?/L U? ` 
The Pe parameter is the fraction of the convection term and the thermal 
diffusion term: 


loou:VT| oocUATL! ocUL | 

|kV?T]| KLAT rk 

The 6 parameter is the ratio of the viscous dissipation term and the convection 
term: 


Pe. 


Veal aU? Wg 
loocu-VT| aoc UAT/L LoocAT 


4.3.3 The Grashof, Prandtl, and Eckert numbers 


The problem with the above dimensionless numbers is that they involve U, 
but U is implied by AT. Assuming that the convection term is much bigger 
than the viscous diffusion term, the momentum equation features a balance 
between the buoyancy term and the convection term: 


|oou:Vul ~ oogß AT. 


Translating this similarity to scales, 


ooU?/L ~ ogb AT, 


gives an U in terms of AT’: 


U=./BLAT. 


The Reynolds number with this U now becomes 
L VgBL3 AT 
Rer = = = gt =Grl/?, 


where Gr is the Grashof number in free thermal convection: 
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The Grashof number replaces the Reynolds number in the scaled equations 
of free thermal convection. We shall soon look at its interpretations, which 
are not as straightforward as for the Reynolds and Peclet numbers. 

The above choice of U in terms of AT results in y equal to unity: 


B gGLAT r gBßLAT = 
Ye ghar 


The Peclet number can also be rewritten as 


L = PrRe * = PrRé;’, 
K K H 


where Pr is the Prandtl number, defined as 


Pr= ie ; 
k 

The Prandtl number is the ratio of the momentum diffusivity (kinematic 
viscosity) and the thermal diffusivity. Actually, more detailed analysis shows 
that Pr reflects the ratio of the thickness of the thermal and velocity boundary 
layers: when Pr = 1, these layers coincide, while Pr « 1 implies that the 
thermal layer is much thicker than the velocity boundary layer, and vice 
versa for Pr œ 1. 

The 6 parameter is in free convection replaced by a combination of the 
Eckert number (Ec) and the Reynolds number. We have that 


and consequently 


Writing 


shows that the Eckert number can be interpreted as the ratio of the kinetic 
energy of the flow and the thermal energy. 

We use Gr instead of Rer in the momentum equations and also instead 
of Pe in the energy equation (recall that Pe = PrRe = PrRer = PrGr~!/?), 
The resulting scaled system becomes 
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Oo pa Va = Vp t Gr 0a Fk + Tk, (4.49) 
V-u=0, (4.50) 
Gri/? (= +ü- v7) = Pr! VPT + 2EcGr 1/74, 58;;. (4.51) 


The Grashof number plays the same role as the Reynolds number in the 
momentum equation in free convection. In particular, it turns out that Gr 
governs the transition between laminar and turbulent flow. For example, the 
transition to turbulence occurs in the range 108 < Gr < 109 for free convection 
from vertical flat plates. Gr is normally interpreted as a dimensionless number 
expressing the ratio of buoyancy forces and viscous forces. 


Interpretations of the Grashof number. Recall that the scaling leading 
to the Grashof number is based on an estimate of U from a balance of the 
convective and the buoyancy terms. When the viscous term dominates over 
convection, we need a different estimate of U, since in this case, the viscous 
force balances the buoyancy force: 


uV?ur~ oogGAT = pU/L? ~ oog SAT, 
This similarity suggests the scale 


B gBL? AT 


V 


U 


Now 


loou-Vul UL gßBAT _ G 


|uV?u] V v = 


The result means that Gr!/? measures the ratio of convection and viscous 
forces when convection dominates, whereas Gr measures this ratio when vis- 
cous forces dominate. 

The product of Gr and Pr is the Rayleigh number, 


Ra= gGL? AT onc 
VK i 
since 
L3 AT L3 AT 
Gir =R pra 2 : ie ae eof = Ra. 
V K VR 


The Rayleigh number is the preferred dimensionless number when studying 
free convection in horizontal layers [2,9]. Otherwise, Gr and Pr are used. 
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4.3.4 Heat transfer at boundaries and the Nusselt and 
Biot numbers 


A common boundary condition, modeling heat transfer to/from the surround- 
ings, is 
OT 


—ky~=A(T'—Ts), (4.52) 


where 0/On means the derivative in the normal direction (n-V), h is an 
experimentally determined heat transfer coefficient, and Ts is the temperature 
of the surroundings. Scaling (4.52) leads to 


-Æ — =h(ATT +T -T;), 


and further to 


OT AL. Ts—To 


T =6(T-T, 
On K ee AT ) ( s), 
where the dimensionless number ô is defined by 
hL 
ô= —, 
k 


and T, is simply the dimensionless surrounding temperature, 


T T's — To 
° AT 
When studying heat transfer in a fluid, with solid surroundings, 6 is known 
as the Nusselt number? Nu. The left-hand side of (4.52) represents in this case 
heat conduction, while the right-hand side models convective heat transfer 
at a boundary. The Nusselt number can then be interpreted as the ratio of 
convective and conductive heat transfer at a solid boundary: 


|n(T—Ts)|_ k 
REIL aD 


The case with heat transfer in a solid with a fluid as surroundings gives the 
same dimensionless 6, but the number is now known as the Biot number’®. It 
describes the ratio of heat loss/gain with the surroundings at the solid body’s 
boundary and conduction inside the body. A small Biot number indicates 
that conduction is a fast process and one can use Newton’s law of cooling 
(see Section 2.1.7) instead of a detailed calculation of the spatio-temporal 
temperature variation in the body. 


=Nu. 


https: //en.wikipedia. org/wiki/Nusselt_number 
Snttps://en.wikipedia. org/wiki/Biot_number 


4.4 Compressible gas dynamics 121 


Heat transfer is a huge engineering field with lots of experimental investiga- 
tions that are summarized by curves relating various dimensionless numbers 
such as Gr, Pr, Nu, and Bi. 


4.4 Compressible gas dynamics 


4.4.1 The Euler equations of gas dynamics 


The fundamental equations for a compressible fluid are based on balance of 
mass, momentum, and energy. When molecular diffusion effects are negligible, 
the PDE system, known as the Euler equations of gas dynamics, can be 
written as 


2 +V-(ou) =0, (4.53) 
Meu) +V- (ouu) =—Vp+ of, (4.54) 
CE LY. (u(H +p) =0, (4.55) 
with E being 
E= pe-+50u-u. (4.56) 


In these equations, u is the fluid velocity, o is the density, p is the pressure, 
E is the total energy per unit volume, composed of the kinetic energy per 
unit volume, Zou -u, and the internal energy per unit volume, ge. 

Assuming the fluid to be an ideal gas implies the following additional 
relations: 


€=CyT, (4.57) 


R 1 
p= oRT= (E— 5ou-u), (4.58) 


Cy 
where cy, is the specific heat capacity at constant volume (for dry air cy = 
717.5Jkg~'K~'), R is the specific ideal gas constant (R = 287.14Jke~'K~'), 
and T is the temperature. 
The common way to solve these equations is to propagate 0, gu, and E 
by an explicit numerical method in time for (4.53)-(4.55), using (4.58) for p. 
We introduce dimensionless independent variables, 
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_. - 2 — y D 13 7 t 
4 Tiel SS Z = — = 
L = y L 3 L a. te >) 
and dimensionless dependent variables, 
soa ee eh ee 
u= 7, OSE T Des E= —. 
U Qc Pe Ee 


Inserting these expressions in the governing equations gives 


ete (ou) =0, 
AGO) 5 Ey (gaa?) =- 1 vps ea, 
a Vv (u(EcE + pep)) = 0, 
p= (BoB ~ S0ctedia ix) 


A natural choice of time scale is te= L/U. A common choice of pressure scale 
is pe = 0-U?. The energy equation simplifies if we choose Ee = pe = o-U?. 
With these scales we get 


ð - 
Sp TV (gu) =0, 
tew) +V (oun!) =-Vp+aðf, 
ðE - z 
oY (u(E+p))=0, 
on Ra i a 
p= —(E- 30u u), 


where a is a dimensionless number: 
_ Ele 
a= Fe 
We realize that the scaled Euler equations look like the ones with dimension, 
apart from the a coefficient. 
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4.4.2 General isentropic flow 


Heat transfer can be neglected in isentropic flow’, and there is hence an 
equation of state involving only @ and p: 


p=F(@). 


The energy equation is now not needed and the Euler equations simplify to 


Oo 7 


Elimination of the pressure. A common equation of state is 


F(0) = po (£). 


where y = 5/3 for air. The first step is to eliminate p in favor of @ so we get 
a system for o and u. To this end, we must calculate Vp: 


y-1 
garoa F(o)=e2 (+) 


where 


is the speed of sound within the fluid in the equilibrium state (see the sub- 
sequent section). Equation (4.60) with eliminated pressure p reads 


a = 
o% +ou-Vut+ee £ Vo=0. (4.61) 
ðt 20 

The governing equations are now (4.59) and (4.61). Space and time are 


scaled in the usual way as 


me x = y = Z 
t=, Yrs 27, i 
L L L te 


The scaled dependent variables are 


Then F’ (o) = cgt. 


Thttps ://en.wikipedia.org/wiki/Isentropic_process 
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Inserting the dimensionless variables in the two governing PDEs leads to 


0c90 , aU ae 
e V -(du) =0 
a Tak ue a 

U ðū oU? _ = 2 
0 + ou Vat E (£) 0" tyo=0. 


The characteristic flow velocity is U so a natural time scale is te = L/U. This 
choice leads to the scaled PDEs 


Oo = 
art (8u) =0, (4.62) 
ðu - iee 
0 tõu Vü+M? (=) o’-!Vo=0, (4.63) 
20 
where the dimensionless number 
vee 
Co 


is known as the Mach number. The boundary conditions specify the character- 
istic velocity U and thereby the Mach number. Observe that (4.63) simplifies 
if 0e = 0o is an appropriate choice. 


4.4.3 The acoustic approximation for sound waves 


Wave nature of isentropic flow with small perturbations. A model 
for sound waves can be based on (4.59) and (4.61), but in this case there 
are small pressure, velocity, and density perturbations from a ground state at 
rest where u = 0, 0 = 00, and p= po = F (00). Introducing the perturbations 
ô= 0-0 and &@, (4.59) and (4.61) take the form 


06 x A 
Ss TV (0+ 0) =0, 


OU RNA A 2 ô oa a 

(00 +ô) + (00+ 0): Va+eG (1+ — Vo=0. 
ot 020 

For small perturbations we can linearize this PDE system by neglecting all 

products of ô and @. Also, 1+ 6/09 ~ 1. This leaves us with the simplified 

system 


4.4 Compressible gas dynamics 125 


06 A 
at V-ú=0 
Ət 20 , 


ot ‘ 
205, + OVO =0. 


Eliminating û by differentiating the first PDE with respect to t and taking 
the divergence of the second PDE gives a standard wave equation for the 
density perturbations: 


Similarly, 6 can be eliminated and one gets a wave equation for &, also with 
wave velocity co. This means that the sound perturbations travel with velocity 
Co. 

Basic scaling for small wave perturbations. Let 0, and uc be charac- 
teristic sizes of the perturbations in density and velocity. The density will 
then vary in [@0 — @c, 00 + 0c]. An appropriate scaling is 


Q2— 20 
Qc 


0 = 
such that g € [—1,1]. Consequently, 


_ & 
20 
Note that the dimensionless @ is expected to be a very small number since 
0c < oo. The velocity, space, and time are scaled as in the previous section. 
Also note that 09 and po are known values, but the scales oe and U are 
not known. Usually these can be estimated from perturbations (i.e., sound 
generation) applied at the boundary. 
Inserting the scaled variables in (4.59) and (4.61) results in 


0= 00+ 0c0=o0(1+a8@), a 


2098 | QU = =~ 
= + V- (14 = 0; 

00 28 WUS (14 aga) 

0U „ôl aU? a- n, 002 ny- 5- 
14 =+ 14 : — co (1 =0. 
i (1+a0) Ai 7 (1+aa)u Vuta co( +a) “Vo=0 


Since we now model sound waves, the relevant time scale is not L/U but 
the time it takes a wave to travel through the domain: te = L/co. This is a 
much smaller time scale than in the previous section because cp >> U (think 
of humans speaking: the sound travels very fast but one cannot feel the 
corresponding very small flow perturbation in the air!). Using te = L/ug we 
get 
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pl +MV-((1+a@)u) =0, 


Ot 
_, OU oes -1 Since 
(1+a0) +M(1+a0)u-Vu+aM*(1+a@)” “Vo=0. 
For small perturbations the linear terms in these equations must balance. 
This points to M and a being of the same order and we may choose a = M 
(implying oe = @9M) to obtain 


Oui g ee 
Sp + Ma Vat (1+Ma)” 293=0. 


Now the Mach number, M, appears in the nonlinear terms only. Letting 
M — 0 we arrive at the following linearized system of PDEs 


oe as 
Ou = 

=a 0 = 4. 
ap + VO=0, (4.65) 


The velocity u can be eliminated by taking the time derivative of (4.64) 
and the divergence of (4.65): 


-= eG, (4.66) 
which is nothing but a standard dimensionless wave equation with unit wave 
velocity. Similarly, we can eliminate @ by taking the divergence of (4.64) and 


the time derivative of (4.65): 


247 = 
a = Vu. (4.67) 


We also observe that there are no physical parameters in the scaled wave 
equations. 


4.5 Water surface waves driven by gravity 


4.5.1 The mathematical model 


Provided the Weber number (see section 4.2.7) is sufficiently small, capillary 
effects may be omitted and water surface waves are governed by gravity. 
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For large Reynolds numbers, viscous effects may also be ignored (except in 
boundary layers close to the bottom or the surface of the fluid). The flow of 
an incompressible homogeneous fluid under these assumptions is governed by 
the Euler equations of motion on the form 


V-u=0, (4.68) 
DY oa ie Gy Pile = 0 (4.69) 
Dt 5 pt+gk=0. ; 
When the free surface position is described as z = n(a,y,t), with z as the 
vertical coordinate, the boundary conditions at the surface read 


o 

S +u-Vn=uw, (4.71) 
where ps is the external pressure applied to the surface. At the bottom, 
z = —h(x,y), there is the no-flux condition 

ðh ry Oh 

az” T w 


In addition to p and g we assume that a typical depth he, a typical wavelength 
Xe, and a typical surface elevation A, which then by definition is a scale 
for 7, are the given parameters. From these we must derive scales for the 
coordinates, the velocity components, and the pressure. 


4.5.2 Scaling 


First, it is instructive to define a typical wave celerity, ce, which must be linked 
to the length and time scale according to Ce = Àc/te. Since there is no other 
given parameter that matches the mass dimension of p, we express ce in terms 
of A, Ac, he, and g. Most of the work on waves in any discipline of physics is 
devoted to linear or weakly nonlinear waves, and the wave celerity must be 
presumed to remain finite as A goes to zero (see, for instance, Section 4.4.3). 
Hence, we may assume that ce must depend on g and either he or Ae. Next, 
the two horizontal directions are equivalent with regard to scaling, implying 
that we have a common velocity scale, U, for u and v, a common length scale 
L for x and y. The obvious choice for L is Ae, while the “vertical quantities” w 
and z have scales W and Z, respectively, which may differ from the horizontal 
counterparts. However, we assume that also the length scale Z remains finite 
as A— 0 and hence is independent of A. This is less obvious for Z than for ce 
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and te, but may eventually be confirmed by the existence of linear solutions 
when solving the equation set. From the linear part of (4.71) and (4.68) we 
obtain two relations between velocity and coordinate scales by demanding 
the non-dimensionalized terms to be of order unity 


A U W 

he? t D AZ 
These relations are indeed useful, but they do not suffice to establish the 
scaling. 

The pressure may be regarded as the sum of a large equilibrium part, bal- 
ancing gravity, and a much smaller dynamic part associated with the presence 
of waves. To make the latter appear in the equations we define the dynamic 
pressure, pg, according to 


(4.72) 


P=Ps— pgz + Pa; 


and the pressure scale pe = pgA for pa then follows directly from the surface 
condition (4.70). 
The equation set will be scaled according to 


L’ L’ gl A’ D” U’ w’? De! 


In the further development of the scaling we focus on two limiting cases, 
namely deep and shallow water. 


4.5.3 Waves in deep water 


Deep water means that he >> Ac. Presumably the waves will not feel the 
bottom, and h as well as he are removed from our equations. The bottom 
boundary condition is replaced by a requirement of vanishing velocity as 
z = —co. Consequently, ce must depend upon Ae and g, leaving us with 
Ce = VgXc and Z = àc = L as the only options. Then, te = \/Ac/g and (4.72) 
implies U = W = cof = eco, where we have introduced the non-dimensional 
number 


= 
which is the wave steepness. The equality of the horizontal and the vertical 
scale is consistent with the common knowledge that the particle orbits in 
deep water surface waves are circular. 


The scaled equations are now expressed with € as sole dimensionless num- 
ber 
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V-a=0, (4.73) 
oe + eu- Vat Vpa=0. (4.74) 


The surface conditions, at z = en, become 


Pa=", (4.75) 


— +u: V= Ù, (4.76) 
while the bottom condition is replaced by 


u—0, (4.77) 


as Z => —oo. 


4.5.4 Long waves in shallow water 


Long waves imply that the wavelength is large compared to the depth: A, >> 
he. In analogy with the reasoning above, we presume that the speed of the 
waves remains finite as A, —> oo. Then, ce must be based on g and he, which 
leads to ce = Vghe and te = Ac /Vghe. The natural choice for the vertical 
length scale is now the depth; Z = he. Application of (4.72) then leads to 
W = cc A/àc and U = c, A/he. 
Introducing the dimensionless numbers 
A he 
he’ HZ> 


C 


Q = 
we rewrite the velocity scales as 


W = ace, U = ace. 


We observe that W « U for shallow water and that particle orbits must be 
elongated in the horizontal direction. 

The equation set is now most transparently written by introducing the 
horizontal velocity üp = ui+ vj and the corresponding horizontal components 
of the gradient operator, Vp: 
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= Ow 
Vai Z =0, (4.78) 
oz 
Ou = Ou = 
CEE aun A Taine a =O: (4.79) 
Ot Oz 
ðr OCH OOD 
2 
Z : 2A = =Q.. 4. 
H E +aup Puto-+ aid) + DE 0 (4.80) 
Surface conditions, at z = an, now become 
Te (4.81) 
On = 
apt th: Va = ü, (4.82) 


while the bottom condition is invariant with respect to the present scaling 


Vp tpn =Ù. (4.83) 


An immediate consequence is that pg remains equal to 7 throughout the water 
column when u? — 0, which implies that the pressure is hydrostatic. The 
above set of equations is a common starting point for perturbation expansions 
in e and u? that lead to shallow water, KdV, and Boussinesq type equations. 


4.6 Two-phase porous media flow 


We consider the flow of two incompressible, immiscible fluids in a porous 
medium with porosity (a). The two fluids are referred to as the wetting® 
and non-wetting fluid. In an oil-water mixture, water is usually the wetting 
fluid. The fraction of the pore volume occupied by the wetting fluid is denoted 
by S(a,t). The non-wetting fluid then occupies 1 —S of the pore volume (or 
(1—S)@ of the total volume). The variable P(æ,t) represents the pressure 
in the non-wetting fluid. It is related to the pressure P, in the non-wetting 
fluid through the capillary pressure pe = Pa — P, which is an empirically 
determined function of S. 

From mass conservation of the two fluids and from Darcy’s law for each 
fluid, one can derive the following system of PDEs and algebraic relations 
that govern the two primary unknowns $ and P: 


Shttps://en.wikipedia.org/wiki/Wetting 
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V-v=—(Qn+Qu); (4.84) 
v=—-A'VP+ Awd, (S)VS + (Aw Ow + Anon) gk, (4.85) 
oF + F(S}U-VS=V-(ha(S)Pl(S)VS)+ 
o fa Qnt Qu) = Qu (4.86) 
Qu = 2, (4.87) 
Ow 
_ In 
Qn = a (4.88) 
(8) = Z kru(3), (4.89) 
An(S) = Ž kra), (4.90) 
MlS) = Aw(S) + An(S), (4.91) 
krul S) = Kum | (4.92) 
b 
krn(S) = Knm eS , (4.93) 
fu(S) = x, (4.94) 
Gw(S) = hw(S) (On — 0w), (4.95) 
hw(S) = —An(S) fw(S). (4.96) 


The permeability of the porous medium is K (usually a tensor, but here 
taken as a scalar for simplicity); 4w and un are the dynamic viscosities of the 
wetting and non-wetting fluid, respectively; Ow and on are the densities of 
the wetting and non-wetting fluid, respectively; qw and qn are the injection 
rates of the wetting and non-wetting fluid through wells, respectively; Swr 
is the irreducible saturation of the wetting fluid (i.e., S > Swr); Snr is the 
corresponding irreducible saturation of the non-wetting fluid (i.e., (1 — S) > 
Snr), Kwn and Ky, are the maximum values of the relative permeabilities 
krw and krn, respectively, and a and b are given (Corey) exponents in the 
expressions for the relative permeabilities. 

The two PDEs are of elliptic and hyperbolic/parabolic nature: (4.84) is 
elliptic since it is the divergence of a vector field, while (4.86) is parabolic 
(hw > 0 because p!.(S) >0 and àn as well as fy are positive since krn > 0 
and krw > 0). Very often, p’, is small so (4.86) is of hyperbolic nature, and 
S features very steep gradients that become shocks in the limit p’, +0 and 
(4.86) is purely hyperbolic. A popular solution technique is based on operator 
splitting at each time level in a numerical scheme: (4.84) is solved with respect 
to P, given S, and (4.86) is solved with respect to S, given P. 
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The saturation S is a non-dimensional quantity, and so are ¢, krw, Krn, 
Kum, Knm, fw, and fl. The quantity v is the total filtration velocity, i.e., 
the sum of the velocities of the wetting and non-wetting fluid. An associated 
velocity scale ve is convenient to define. It is also convenient to introduce 
dimensionless fractions of wetting and non-wetting fluid properties: 


0 = 0w, 
On 
On = QQ, a=—_, 
Ow 
H = Hw, 
H 
Un = HP, =, 
Hw 
q 
Q= w=, 
Hi q 
n=Q-, YS a 
a dw 


AlS) = Awl S) AnS) = Acà, At = àw te lAn. 


As we see, Àc is the characteristic size of any “lambda” quantity, and a bar 
indicates as always a dimensionless variable. The above formulas imply 


hw(S) = —AcB~*An(S) fuls),  Gw(S) = hw(S)o(a—1). 
Furthermore, we introduce dimensionless quantities by 


_ £ _ v P P _ De 
4 = — v = — = — = —. 
T ve’ po PeT p 


Inserting the above scaled quantities in the governing PDEs results in 
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v =- ata), (4.97) 
v= Pore x VP4 Aches 5 (8) 7S+ 

veL Ve 

DÀ Ay tap Än)k, (4.98) 
Os teVe ! a _ tePere —1 / 
bag + fel VS = SIV: OE AlS) fuwl) PSV S) 

teg OGw aa 

Be ttefuQll+aT")-teQ. (4.99) 


As usual, L is taken as the characteristic length of the spatial domain. Since 
Ue is a velocity scale, a natural time scale is the time it takes to transport a 
signal with velocity ve through the domain: te = L/ve. The diffusion term in 
the equation (4.102) then gets a dimensionless fraction 


LP.Xe 
UcL2 ` 


Forcing this fraction to be unity gives 


P, 
Ye= Ao. 


We realize that this is indeed a natural velocity scale if the velocity is given 
by the pressure term in Darcy’s law. This term is K/j times the pressure 
gradient: 


KP. 


K 
—|VP| ~ 
H u L 


= ewe = Ve. 
We have here dropped the impact of the relative permeabilities Aw or An 
since these are quantities that are less than or equal to unity. 

The other term in Darcy’s law is the gravity term that goes like Acog 
(again dropping relative permeabilities). The ratio of the gravity term and 
the pressure gradient term in Darcy’s law is an interesting dimensionless 
number: 


_ Aceg_ _ Log 
NePe/E Pe 
This number naturally arises when we discuss the term 
teg 0G nas! mT. wi xt 1 Os 
of TER = (a — 1)8 TNCS) fur(S) + An lS) AlS) S 


Introducing another dimensionless variable, 


L? 
c=teQ= re, 
cle 
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we can write (4.97)-(4.99) in the final dimensionless form as 


V-0=—-e€(1+a+y), (4.100) 
B= uV P+ Awd (S)VS +5 Aw tab An)k, (4.101) 
oS + $1,(5)B- VS = -9 (8n (8) ful 5) PCS?) 
(a= 1)8-14(%% (5) fu(S) + 7n () A) 
eful +a™ty)—e. (4.102) 


The eight input parameters L, qw, qn, Hw, Hn, Ow; On, and K are reduced 
to five dimensionless parameters a, 8, y, 6, and e. There are six remaining 
dimensionless numbers to be set: Kwm, Knm, Swr, Snr, a, and b. 
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